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FOREWORD 


This final report describing the formulation of the Six- 
Degree-of -Freedom Program to Optimise Simulated Trajectories 
(6D POST) is provided in accordance with Part 3.0 of NASA Con- 
tract NAS 1-13300. The report is presented in three volumes as 
follows: 

Volume I - 6D POST - Formulation Manual; NASA CR- 132741 

Volume II - 6D POST - Utilisation Manual; NASA CR— 1 32742 

Volume III - 6D POST - Programmer's Manual. NASA CR-132743 

This work was conducted under the direction of Mr. Howard Stone 
and Mr. Richard Powell of the Space Systems Division, National 
Aeronautics and Space Administration, Langley Research Center. 
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SUMMARY 


This report documents the basic equations and models used 
In the slx-degree-of-freedom version of the program to optimize 
simulated trajectories (CD POST). 

6D POST, a direct extension of the point mass version of 
POST, Is a general purpose rigid body slx-degree-of-freedom pro- 
gram. The program can be used to solve a wide variety of 
atmospheric flight mechanics and orbital transfer problems for 
powered or unpowered vehicles operating near a rotating oblate 
planet. The principal features of 6D POST are: an easy to use 
NAMELIST type input procedure, an Integrated set of Flight Con- 
trol System (FCS) modules, and a general-purpose discrete 
parameter targeting and optimization capability. 

6D POST Is written In FORTRAN IV for the CDC 6000 series 
computers . 

Other volumes In the final report are: 

Volume II - Utilization Manual - Documents Information 
pertinent to users of the program. It describes the 
Input required and output available for each of the 
trajectory and targetlng/optlmlzatton options. 

Volume III - Programmers Manual - Documents the program 
structure and logic, subroutine descriptions, and other 
pertl nent programml ng 1 nformatl on . 
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I. INTRODUCTION 


The six (6 )-degree-of- freedom program to optlnize simulated 
trajectories is a general purpose FORTRAN program for simulating 
rigid body trajectories of aerospace type vehicles. The program 
can be used to solve a wide variety of performance, guidance, and 
flight control problems for atmospheric and orbital vehicles. For 
example, typical applications of 6D POST Include: 

1) Guidance and flight control system simulation and analysis; 

2) Loads and dispersion type analyses; 

3) . General-purpose 6D simulation of controlled and uncon- 

trolled vehicles; 

4) 6D performance validation. 

One of the key features of 6D POST Is an easy to use NAMELIST- 
type input procedure. This feature significantly reduces input 
deck set-up time (and costs) for 6D studies that require the normal 
large amount .of input. data. In addition., the general applicability 
of 6D POST is further enhanced by a general-purpose discrete param- 
eter targeting and optimization capability. This capability can 
be used to solve a- broad spectrum of problems related to the Impact 
of the control system design on the performance characteristics of 
aerospace vehicles. 

The „ ic simulation flexibility is achieved by decomposing 
the trajectory into a logical sequence of simulation segments. 

These trajectory segments, referred to as phases, enable the tra- 
jectory analyst to model both the physical and the nonphysical 
aspects of the simulation accurately ani efficiently. By segment- 
ing the mission into phases, each phase can be modeled and simulated 
in a manner most appropriate to that particular flight regime. For 
example, the planet model, the vehicle model, and the simulation 
options can be changed in any phase to be compatible with the level 
of detail required in that phase. 

Every computational routine in the program can be categorised 
according to five basic functional elements. These elements are: 
the planet model, the vehicle model, the trajectory simulation 
model, the auxiliary calculations module, and the targeting and 
optimisation module. The planet model is composed of an oblate 
spheroid model, a gravitational model, an atmosphere model, and 
a winds model. These models define the environment in which the 
vehicle operates. The vehicle model comprises mass properties, 
propulsion, aerodynamics and aeroheating, an airframe model, a 
navigation and guidance model, and a flight control sy^cem model. 
These models define the basic vehicle simulation characteristics. 

The trajectory simulation models are the event-sequencing module 


that controls the program cycling, table interpolation routines, 
and several standard numerical Integration techniques. These 
models arf, used in numerically solving the translational and ro- 
tational equations of motion. The auxiliary calculations module 
provides for a wide variety of output calculations. For example , 
conic parameters, range calculations, and tracking data are among 
the many output variables computed. The targeting- and optimiza- 
tion module provides a general discrete parameter iteration capa- 
bility. The user can select the optimization variable, the de- 
pendent variables, and the Independent variables from a list of 
more than 400 program variables. An accelerated projected 
gradient algorithm is used as the basic optimization technique. 
This algorithm is a combination of Rosen's projection method for 
nonlinear programming and Davidon's variable metric method for 
unconstrained optimization. In the targeting mode, the minimum 
norm algorithm is used to satisfy the trajectory constraints. 

The cost and constraint gradients required by these algorithms 
are computed as first differences calculated from perturbed tra- 
jectories. To reduce the costs of calculating numerical sensi- 
tivities, only that portion of the trajectory influenced by any 
particular independent variable is reintegrated on the perturbed 
runs. This feature saves a significant amount of computer time 
when targeting and optimization is performed. 

Basic program macroiogic is outlined in figure I— 1, which 
illustrates the linkage between the simulation and the iteration 
modules . 
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II. LIST OF SYMBOLS AND JkBBRBKIATIPN.S. 


Math symbol 

a 

-AB " ( A AXB» A AYB» A AZb) 
[AB] 

V 

*h* a mp» 

M 

*s 

^SB " ( a sxb* a syb* a szb) 

^SI “ (^XI' ^YI’ ^Zl) 

^TB “ ( A I»' ^TYB* A TZb) 

^L 

*21* k ZBL* A ZA 


Internal Fortran symbol 
SEMJAX 

AB(I) 

AR 

AMXB, AMYB, AMZB 

A(I) 

ASM 

AXB, AYB, AZB 
AS XI, ASYI, ASZI 

AZL 

AZVELI, AZVELR, 
AZVBLA 


Definition 

semimajor axis, m (ft) 

aerodynamic acceleration 
In the body frame, mps 2 
(fps 2 ) 

matrix transformation from 
the A-frame to the B-frame 

nosale exit area of each 
rocket engine, m 2 (ft 2 ) 

constants 

total aerodynamic moment 
about the roll, pitch, yarn 
axes, N-m (ft-lb) 

Davldon deflection matrix 
component 

total sensed acceleration, 
mps 2 (fps 2 ) 

total- sensed acceleration 
in the body frame, mps 2 
(fps 2 ) 

total sensed acceleration 
In the Inertial frame, 
mps 2 (fps 2 ) 

thrust acceleration In the 
body frame, mps 2 (fps 2 ) 

azimuth of the a, axis, 
rad (deg) 

aslmuth of the Inertial, 
relative, and atmospheric 
relative velocity vectors, 
rad (deg) 


II-l 




~T 


V 


AZWT 
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wind azimuth , rad (deg) 

azimuth of the slant range 
vector to the tracking 
station, rad (deg) 

boundary for 1 con- 
straint 

Davidon deflection matrix 

boundary of region R 

local boundary hyper- 
surface 

axial, side force, and 
normal aerodynamic force 
coefficients 

component of C A* V S 

that is not multiplied 
by a mnemonic multiplier 

drag and lift coefficients 

drag and lift coefficient 
components that are not 
multiplied by a mnemonic 
variable 



Math symbol 

Internal Fortran symbol 

Definition 

V C n 

CM, CW 

pitch and yaw moment co- 
efficients 

C S 

CS 

speed of sound, mps (fps) 

C(u) 

E(I) 

constraint functions 

D- 

DRAG 

aerodynamic drag, M (lb) 

E 

ECCAN 

eccentric anomaly 

[E] 


Euler parameter matrix 

JB 

ECCEN 

eccentricity 

- “ (V V * 2 * *3) 

E0(I) 

Euler parameters 

II 

E(I) 

active constraint error 
vector 

* 

£ 

WE(I) 

weighted error vector 

P 

0PTVAR 

optimisation function 

£ 

— 

nonlinear vector-valued 
function 

^AB " ( F AXB’ F AYB* F AZB 

J FAXB, FAYB, FAZB 

aerodynamic forces In the 
body frame, N (lb) 

^tb “ ( f txb* f tyb* f tzb) 

| FTXB, FTYB, FTZB 

thrust forces In the body 
frame, N (lb) 

[GA] 

GA(I) 

matrix transformation from 
the G- frame to the A-frame 

£t - (°xr Sr Si) 

GX1, GYI, GZ1 

total gravitational 
acceleration In the EC1- 
frame, mps 2 (fps 2 ) 


Mcth symbol 


Internal Fortran symbo l Definition 



DG(l) 

difference In the gradient 
vector VF between the cur- 
rent and previous Itera- 
tion 

H 

— 

gravitations const, nt 

h 

ALTIT0 

oblate altitude. .■ 

h - (h XI , hyj, 

ANGM0M 

*n\‘ nomtei. . ups 2 

V h P 

ALTA, ALTP 

altitude of apogee and 
perigee, km (n mi) 

h b 

HB 

base altitude used in 
atmospheric calcul at ions , 
m (ft) 

h 

c 

P2 

constraint function 

H g 

HT 

geopotential altitude, m 
(ft) 

h e 

R i 

— 

heating ratios 

h T 

TRKHT1 

altitude of tracker, m 
(ft) 

*o 

P1NET 

estimated net cost func- 
tion 

1 

INC 

relative-frame orbital 
inclination, rad (deg) 

UB] 

IB (I) 

matrix transformation from 
the ECl-frame to the body 
frame 

UGl 

IG(I) 

matrix transformation from 
the ECI-frame to the geo- 
graphic frame 
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Internal Fortran symbol. Definition 


Math symbol 


[IL] 

IL(I) 

[IP] 

IP (I) 


ISPV 

I «P 


[J] 

AC0B(J) 

J 2» J 3* J 4 

J2, J3, J4 

k ■ ^2* ^3 * ^ 4 ) 

... 

K i 

— 

L 

LIFT 

[LB] 

LB(I) 

ft 

LREF 

M 

MACH 

M 

MEAN 

M 

— 

m 

MASS 

M f 

— - 

m 

NAC 

n a 



matrix transformation from 
the ECI-frame to the 
launch frame 

matrix transformation from 
the ECI-frame to the 
planet frame 

rocket specific impulse , 
s 

constraint Jacobian matrix 
gravitational constants 

Runge-Kutta constants 


constants 

aerodynamic lift* N (lb) 

matrix transformation 
from the launch frame to 
the body frame 

aerodynamic reference 
length, m (ft) 

Mach number 

mean anomaly, rad (deg) 

pitch and yaw moment equa- 
tions 

vehicle mass, kg (slug) 

mnemonic table multiplier 
for table f 

number of active con- 
straints 
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Math symbol 

Internal Fortran symbol 

Definition 

n c 

NDEPV 

number of constraints 

[PI. ih 

PR0J (I) 

projection operators used 
in the projected gradient 
method 

P(h) 

PRES 

atmospheric pressure, 
N/m 2 (psf) 

P 1 

PI 

weighted optimization 
variable 

P 2 

P2 

weighted constraint error 
function 

Q 

TLHEAT 

total heat, J/m 2 (Btu/ft 2 ) 

"<r 

DYNP 

dynamic pressure, N/m 2 
(lb/ft 2 > 

• ♦ 
^lam* ^turb 

HEATRT, HTURB 

laminar and turbulent heat 
rate, W/m 2 /s (Btu/ft 2 /s) 

Q(«) » §(u) 

— 

linear manifold and its 
orthogonal . complement 

r a 

RTASC 

right ascension of out- 
going asymptote, rad (deg) 

r a 

AP0RAD 

apogee radius, m (ft) 

[RB] 

RB(I) 

matrix transformation from 
the body reference frame 
to the body frame 

*D 

DPRNG1 

dot-product range, km 
(n mi) 

V *p 

RE, RP 

equatorial and polar 
radius, m (ft) 

R.. 

RN 

nose radius, m (ft) 


Math symbol 


Internal-Fortran symbol Definition 


*NU 

REYN0 

Reynolds number 

- {«j. y v «j) 

XI, XI, 21 

Inertial radius vector 
from center of planet *o 
the vehicle, m (ft) 

r I 

GCRAD 

geocentric radius, m (ft) 

r p 

PGERAD 

perigee radius, m (ft) 

*a 

RS 

radios to oblate surface, 
m (ft) 

Ssr 


slant range vector, m (ft) 

SsRB 


slant range vector In geo- 
graphic frame, m (ft) 

•Sir 

— 

radius vector to tracking 
station, m (ft) 

£ 

S(I) 

direction of search 

s c 

— 

direction of search to 
satisfy the constraints 

e 

iOSS^ 

SL0SIJ 

space losses for tracking 
stations, dB 

S ref 

SSEF 

aerodynamic reference 
area, m 2 (ft 2 ) 

a 0 

— 

direction of* search for 
optimization 

X 

ATEM 

atmospheric temperature, 
•K (°F) 

t 

TIME 

time, s 

V 

— 

jet engine thrust, N (lb) 
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Math symbol 


Internal Fortran symbol Definition 


T n (y) 

••••• 

denotes n fc ^" order table 
Interpolation on the 
variable y 

t r 

THRUST 

total rocket thrust for 
all engines, N (lb) 

t r 

R i 


total resultant rocket 
thrust for engine 1, N 
(lb) 

T 

vac 

TVAC 

vacuum thrust for rocket 
engines, N (lb) 

T 

SP 

TIMSP 

time since perigee pas- 
sage-, s 

T 

TP 

TIMTP 

time to next perigee pas- 
sage, s 

H 

- — 

gravitational potential 
function 

u 

U(I) 

independent variable 
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- (V V “,) 


OB, VB, UB 


components of the atmos- 
pheric relative velocity 
vector expressed in the 
body frame, mps (fps) 





— 

Au 

DU(I) 

v . 

APVEL 

ho 

UA, VA, WA 

hi 

VAXI, VAYI, VAZI 

h m ( V XI» V Yt» V 2l) 

VAI, VYI, VZI 


VELI 

“IC- 

Uj v, W 

h 

VELR 

ho 

UR, VR, MR 

hi “ ( V m» V RYI* V R2l) 

VRXI, VRYI, VR2I 


unit vector along. the. 
radius vector 

unit vector along the 
velocity vector 

change in the independent 
variables 

inertial velocity at 
apogee, mps (fps) 

atmospheric relative 
velocity in the 6-frame, 
mps (fps) 

atmospheric relative 
velocity vector in the 
inertial frame, mps (fps) 

inertial velocity vector 
and its magnitude, mps 
(fps) 

magnitude of V^., mps (fps) 

Inertial velocity in the 
G- frame, ops (fps) 

relative velocity, mps 
(fpe) 

relative velocity in the 
G-frame, mps (fps) 

relative velocity vector 
in the inertial frame, 
mps (fps) 






1 


* 

I 

) 

I 


Math symbol 

Internal Fortran symbol 

Definition 

hi ■ (W W • v w*i) WXI - wzl 

0 

wind velocity vector in 
the inertial frame, tnps 
(Cps) 

\ 

VW 

wind velocity, mps (fpa) 

^WG 

UW, VW, WW 

wind velocity vector in 
the G-frame, mps (fps) 

V 

P 

PGVEL 

perigee velocity, ops 
(fps) 

V 

HYPVEL 

outgoing asymptote 
velocity, mps (fps) 

* 

w 

WD0T 

total time rate of change 
of vehicle weight, N/s 
(lb/s) 

w c 

WEIC0N 

total weight of propel- 
lant consumed, N (lb) 

W G 

• WEIGHT 

gross vehicle weight, N 
(lb) 

W jett 

WJETTM 

jettison weight, N (lb) 

W PC 

— 

weight of propellant con- 
sumed per phase, N (lb) 

\ 

— 

initial propellant weight, 
N (lb) 

w®** 

W — 

maximum flowrate for the 
i th engine, N/e (lb/s) 

W PR 

WPR0P 

weight of propellant re- 
maining, N (lb) 

W 8tg 

WGTSG 

vehicle stage weight, N 
(lb) 

N' N* W 

WU, W0PT, WE 

weighting matrices for u, 
and e 
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H 


Math symbol 

Internal Fortran symbol 

Definition 

*B» y B* *B 

— 

coordinate axes of the 
body frame 

X BR* y BR* Z BR 

— 

coordinate axes of the 
body reference frame 

X cg* y cg’ 2 cg 

XCG, YCG, ZCG 

coordinates of the center 
of gravity in the body 
reference system, m (ft) 

V V *G 


components of a vector in 
the geographic frame, m 
(ft) 

*!» y x » 

XI, YI, ZI 

components of the radius 
vector in the inertial 
frame, m (ft) 

x i 

— 

general state variable 

V y L* *L 

- — 

coordinate axes of the 
launch frame 


. GINTJ 

state vector at the n^ 
event 

*R* y R’ Z R 

— 

components of the radius 
vector in the planet 
frame, m (ft) 

x ref ' y ref* *ref 

XREF, YREF, ZREF 

coordinates of the aero- 
dynamic reference point 
in the body reference 
system, m (ft) 

2. 

DGENV 

general dependent variable 

®* 0 » ® 

ALPHA, BETA, 
BHKAMG 

aerodynamic angle of at- 
tack, sideslip, and bank, 
rad (deg) 

®T 

ALPT0T 

total angle of attack, 
rad (deg) 
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GAMMA! , GAMMAR, 
GAMMAA 

inertial, relative, and 
atmospheric relative 
flight path angles, rad 
(deg) 

GAMMA(I) 

step-size parameter on the 
j th trial step 

— 

increment in eccentric 
anomaly, rad (deg) 

— 

increment in altitude, m 
(ft) 

DT 

increment in time or inte- 
gration step size, s 

DV 

increment in velocity, 
mps (fps) 

VIDEAL 

ideal velocity, mps (fps) 

DLR 

atmospheric velocity loss, 
mps (fps) 

DVCIR 

velocity required to 
circularize an orbit, 
mps (fps) 

DVEXS 

excess velocity, mps (fps) 

GLR 

gravity loss, mps (fps) 

DVMAR 

velocity margin, ups (fps) 

ATLR 

atmospheric pressure 
loss, mps (fps) 

TVLR 

thrust vector velocity 
loss, mps (fps) 


RTASC 


right ascension, rad (deg) 


Math symbol 
n 

Internal Fortran symbol 
ETA 

Definition 

engine throttling param- 
eter 

e 

L0NG 

planet relative longitude 
rad (deg) 

8* 

— 

longitude reference, rad 
(deg) 

9 I 

L0NGI 

inertial longitude, rad. 
(deg) 

V V *ZL 

L0NL, LATL, AZL 

longitude, latitude, and 
aslmuth of t-frame, rad 
(deg) 

^max 

TRUNMX 

maximum true anomaly for 
hyperbolic orbit, rad 
(deg) 

0 T 

X l 

TRKLNI 

longitude of tracker 1, 
rad (deg) 

X 

A2REF 

azimuth reference, rad 
(deg) 

X 

STPMAX 

maximum admissible step 
size for the iteration 
algorithm 

M 

MU 

gravitational constant, 
m 3 /s 2 (ft 3 /s 2 ) 

V 

___ 

index 


11-13 


Math symbol 


Internal Fortran symbol 


GCLAT 


GDLAT 


’i' V 6 I 


R0LI, YAWI, PIT1 


^R* 9 R» *R 


YAWR, PITR, R0LR 


0MEGA 


v V W *J 


R0LBD, PITBD, 
YAWBD 


R0LBDO , PITBDD, 
YAWBDD 
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argument of vehicle (l*e.»- 
angular location of ve- 
hicle, measured from 
ascending node in orbital 
plane), rad (deg) 

atmospheric density, kg/m 3 
( slug/ ft 3 -)- • 

trajectory propagation 

geocentric latitude, rad 
(deg) 

geodetic latitude, rad 
(deg) 

inertial roll, yaw, and 
pitch measured as positive 
rotations from the L-frame, 
rad (deg) 

relative yaw, pitch, and 
roll, measured in a posi- 
tive sense from the geo- 
graphic frame, rad (deg) 

longitude of ascending 
node, rad (deg) 

angular rotation rate of 
planet about the polar 
axis, rad/s (deg/s) 

argument of perigee, rad 
(deg). - 

inertial angular velocity 
components about the body 
axis, rad/s (deg/s) 

inertial angular acceleration 
components about the body 
axis, rad/s 2 (deg/s 2 ) 


Math symbol 


Internal Fortran symbol Definition 


<>A 

refers to atmosphere rela- 
tive variables 

<>C 

refers to center of 
gravity 

< >1 

refers to inertial 
variables 

<>» 

refers to event 

<>, 

refers to thrust applica- 
tion 

<>* 

refers to Earth-relative 
variables 

<>«.£ 

refers to aerodynamic ref- 
erence point 

( ^SL 

refers to sea-level condi- 
tions 

<>vec 

refers to vacuum condi- 
tions 

<>W 

refers to wind relative 
variables 

( )* . 

refers to state from which 
dovnrange and crossrange 
are referenced; refers to 
optimal conditions 

(J 

denotes vector quantity 

( )' 

denotes transpose of a 
vector 

<‘> 

denotes total derivative 
with respect to time 
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Math symbol 

( > + 

( )" 


e 

n 


u 

(A:B) 

3 


I nternal Fortran symbol Dei inltion 

denotes occurrence at the 
positive side of an event 

denotes occurrence at the 
negative side of an event 

is a member of 

intersection of 

union of 

set 

such that 


0 


addition operator 


The following symbols have been added to the six-degree-of -freedom 
portion of the program. 


Math symbol 


Internal Fortran symbol Definition 


C Afia* C A6e* C AfiR* 
C A6f t 


C Dfia* C D6e* C DfiR* 
C D6f 1 


C £6a* C f fie* C ffiR, 
C ifif 1 


CAD A, CADE, CADR, 
CAF(I) , 1-1,3 


CODA, CODE, CDDR, 
CDF (I) , 1-1,3 


CLLDA, CLLDE, CLLDR 
CLLF(l) , 1-1,3 


Incremental axial force 
coefficient per rad (deg) 
for the aileron, elevator, 
rudder, and general 
def lector 

Incremental drag force 
coefficient per rad (deg) 
for the aileron, elevator, 
rudder, and general 
deflector 

Incremental rolling 
moment force coef- 
ficient per rad (deg) 
for the aileron, elevator, 
rudder, and general 
del lector 


Math symb ol 


Intgrnal Fortran ay mbol Definit ion 


C Ua* S.’-a* C L 
c ui 4 

CLUA, CLDt' , CUIR, 
CLF(l), 1-1,3 

mc.j’ Sn.R* 

c m\f . 

1 

CMDA, CMDE, OMDR, 
CMF(l), 1-1,3 

C ni'a* C nl-o* C n,U‘ 
C n6f 1 

CWDA, CWDS, CUT)R, 
CWF(l), 1-1,3 

Slda* Sl6e* Sl5R* 
C M6f 1 

CNDA, CNDE, CNDR, 
CNF (I) , 1-1,3 

Sffia* Sfde* Sfdr’ 

CYDA, CYDE, CYDR, 
CYF(l) , 1-1,3 

d R* V S 

DREFR, DREFP, 
DREFY 

h 

FTTXB(I) , 1-1,3 

hcB 

RCSFXB(I), 1-1,3 

h* 

FTXB(l), 1-1,3 


Increment .*1 lift tore© 
coefficient pc-r rad (deg) 
for the aileron, elevator, 
rudder, and general 
deflector 

Incremental pitching 
moment force coef- 
ficient per rad (deg) 
for the aileron, elevator, 
rudder, and general 
deflector 

Incremental vnvinp moment 
force .uef fi<*U;nt per 
rad (deg) for the aileron, 
elevator, rudder, and 
genera) deflector 

Incremental normal force 
coefficient per rad (deg) 
for the aileron, elevator, 
rudder, and general 
deflector 

Incremental side force 
coefficient per rad (deg) 
for the aileron, elevator, 
rudder, and general 
deflector 

Reference lengths for roll 
pitch, and yaw aerodynamic 
moment coefficients m (ft) 

Total of all nongravlta- 
tlonal forces calculated 
in the body frame, N (lb) 

Force of the RCS engines 
in the body frame, N (lb) 

Total force due to the 
non-RCS engines calculated 
in the body frame, N (lb) 


I 


\ 


Hath symbol 

internal Fortran symbol 

Definition 

Six’ 1 yy* *zz 

IXX, IY V , IZZ 

Momenta of Inertia 
about the body axia 
system 

*XX» *YY* *ZZ* 
*XY* ^YZ* *XZ 

1 XXD , IYYD, IZZD, 
IXYD, IYZD, IXZD 

Time derivations of the 
moment ~ and products 
of Inertia 

I XY’ I xz* *YZ 

IXY, IXZ, IYZ 

Products of Inertia 
about the body axis system 

[K] 

KRDP(l), KPDP(I), 
KRDY(I), KYDY(I) 

Roll nozzle deflection 
matrix 

S'a* K Poe* V'r 

KPDA, KPDE, KPDR 

Mixing logic gains for the 
aileron, elevator, and 
rudder about the y-body 
(pitch) axis 

K R6a* K R6e’ K R*>r 

r‘KDA, KRDE, KRDR 

Mixing luglc gains for 
the aileron, elevator, 
and rudder about the 
x-body (toll) axis 

Sfia’ *Sf6«* S.sr 

KYDA, KYDE, KYDR 

Fixing logic gains for 
the z-body (yaw) axis 

IM] 

— 

Matrix representation 
of the mixing gains 

-AB 

AtttB(i), 1*1,3 

Externa:, moment due to 
the aerodynamic forces, 
N-m (ft-lb) 

55b 

TTMXB(I), 1-1,3 

Total external moment- 
due to thrust, RCS, and 
aerodynamic forces, 

N-m (ft-lb) 

55rcs 

RCSMXBU), 1-1,3 

Net moment due to the RCS 
forces, N-m (ft-lb) 

^TB 

TMXB(l), 1-1,3 

External moment due to 
thrust forces, N-m (ft-lb) 
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Math symbol 

P'» q'. r - 


6.* £o 


6a, 6e, 6 r , 6 f 




i 


60 _ « ( 5 60 , 6 $) 


♦ 

^^AB 




1 


8P* 


'gp' *8P 


Internal Fortran symbol 
PND, QND, RND 

ALP Hi , BETAI , BANKI 


DELA, DELE, DELR 
DELF(I), 1-1,3 

DER(l), DEY(I) 

YAWAC , PITAC, R0LLAC 
DXR(I), 1-1,3 

DXP(I) , DYP(I) , DZP(I) 
GXP, GYP, GZP 


Definition 

Nondimens Iona 1 roll, 
pitch, and yaw body 
rates 

Attitude reference 
angles measured in the 
same sense and order as 
the aerodynamic angles 
but with respect to the 
inertial velocity vector, 
rad (deg) 

General vector represent- 
ing an engine or aero- 
dynamic control surface 
deflection, r«*a (deg), 
the subscript denotes 
the null value. 

Deflection angles for the 
aileron, elevator, rudder, 
and general aerodyuamlc 
control surfaces rad (deg) 

Pitch and yaw glmbal angles 
for the 1-th engine, rad 
(deg) 

Yaw, pitch, and roll auto- 
pilot commands, rad (deg) 

Vector difference between 
the center of gravity 
and the reference point 
for the aerodynamic forces 
(usually the aerodynamic 
center of pressure) , m (ft) 

Vector difference between- 
the center of gravity and 
the engine gimball point 
for the i-th engine, m (ft) 

Location of engine gimbal 
in body reference system 
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III. COORDINATE systems 


6D POST uses numerous coordinate systems to provide the neces- 
sary reference systems for calculating required and optional data. 
These coordinate systems and the key tranaformations are described 
below. 


Coordinate System Definitions 

Earth-centered inertial (ECI) axes (xj, y x , *” This 8y8 “ 
tem is an Earth-centered Cartesian system with z x coincident 
with the North Pole, x x coincident with the Greenwich Meridian 
at time zero and in the equatorial plane, and y x completing a 

right-hand system. The translational equations of motion are 
solved in this system (fig. XXl-1)* 

Earth-centered rota ting (ECR) axes |x^, y R > 2 r) *~ This sys- 
tem is similar to the ECI system except that it rotates with the 
Earth so that x_ is always coincident with the Greenwich Merid 

ian (fig. III-l) . 

Earth position coordinates (<j> g » 6 > ♦- These are the fa- 

miliar latitude* longitude, and altitude designators . Latitude 
is positive in the Norther; Hemisphere. Longitude is measured 
positive East of Greenwich. Altitude is measured positive above 
the surface of the planet (fig. III-l). 

Geographic (G) axes (xg. y G , s Q ). - This system is located 

at the surface of the planet at the vehicle’s current geocentric 
latitude and longitude. The x Q axis is in the local horizontal 

plane and points North, the y G axis is in the local horizontal 
plane and points East, and z G completes a right-hand system. 
This system is used to calculate parameters associated with azi- 
muth and elevation angles (fig. III-2). 

Inertial launch (L) axes (x^ y L » « L ) » “ Thl ® i8 an iner " 

tial Cartesian system that is used as an inertial *®*®'®J?® 

system from which the inertial 4 .® tt X a U au t onmtically located at the 
measured. This coordinate system is automatically 1 





geodetic latitude and Inertial longitude of the vehicle at the 
beginning of the simulation unless overridden by user input of 
LATL and L0NL. The azimuth, A^, Is zero unless overldden by 

user Input. The orientation of this system is such that x^ is 
along the positive radius vector If Is input as the geocen- 

tric latitude, or along the local vertical if 0^ is not input 
or is input as the geodetic latitude, z^ is in the local hori- 
zontal plane and is directed along the azimuth specified by A^, 
and y. completes a right-hand system. This system is intended 

for use in simulating ascent problems for launch vehicles that 
use either inertial platform or strapdown-type angular commands. 
The inertial angles, v are always measured with 

respect to this system and are automatically computed regardless 
of the steering option (IGUID) being used (fig, III-2), 

Body (B) axes (xg, y B , Zg).- The body axes form a right- 

hand Cartesian system aligned with the axes of the vehicle and 
centered at the vehicle’s center of gravity. The Xg axis is 

directed forward along the longitudinal axis of the vehicle ,. yg 
points right (out the right wing), and Zg points downward, com- 
pleting a right-hand system.. All aerodynamic and thrust forces 
are calculated in the body system. These forces are then trans- 
formed to the inertial (I) system where they are combined with 
the gravitational forces (fig, III-3) 



Figure III-3.- Body Frames 


Body reference (BR) axes (x BR , y BR » z BR j The body reference 

system Is a right-hand Cartesian system aligned with the body axes 
as follows. The Xg R axis Is directed along the negative x B axis, 

the y BR axis Is directed along the positive y B axis, and the z BR 

Is directed along the negative z B axis. This system Is used to 

locate the vehicle's center of gravity, aerodynamic reference point, 
and engine glnfcal locations for the static trim operation (fig. 

I 11-3) . 


Orbital elements (h.., h n , 1, n, e, u>).- This Is a nonrectangular 

coordinate system- used In describing orbital motion. The orbital 
elements are apogee altitude, perigee altitude. Inclination, longi- 
tude of the ascending node, true anomaly, and argument of perigee. 

The apogee and perigee altitudes replace the standard orbital ele- 
ments of semimajor axis and eccentricity (ftg. I I 1-4). 



Figure III-4.- Orbital Parameters 



Attitude Angles 


The program contains the following standard attitude refer- 
ence systems: 


1) Inertial Euler angles; 

2) Relative Euler angles; 

3) Aerodynamic angles; 

4) Inertial aerodynamic angles; 

These variables are defined and Illustrated below: 


1 ) 


♦l 


Inertial Euler angles (fig. III-5): 

- Inertial roll angle. The roll 

angle with respect to the L- 
frame (first rotation), y 

la 

- Inertial yaw angle. The yaw 
angle with respect to the L- 
frame (second rotation) , 



y B 


*L 


*B 


0_ - Inertial pitch angle. The pitch 
1 angle with respect to the L- Figure III-5,- 
frame (third rotation) ; 


Inertial Euler Angles 


2) Relative Euler angles (fig. III-6): 


4 > 


R 


♦r 


Relative yaw angle. This is 
the azimuth angle of the 



axis measured clockwise from 
the reference direction (first 
rotation) , 

Relative pitch angle. This is 
the elevation angle of the Xg 

axis above the local horizontal 
plane (second rotation) , 



Relative roll angle. This is 

the roll angle about the Xg Figure III-6.- Relative Euler Anglee 


axis (third rotation) . 
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Transformations 


Numerous matrix transformations are required to transform 
data between the coordinate systems described in the previous 
section. The most Important of these transformations is the [IB] 
matrix. The inverse ^transpose) of this matrix is used to trans- 
form accelerations in the body frame to the planet-centered In- 
ertial frame. The remaining transformations are generally used 
to either compute [IB] or to transform auxiliary data into some 
convenient output coordinate system. 

the [IB] matrix is functionally dependent on the attitude 
of the vehicle. This dependence is described by equations re- 
lated to the attitude steering option selected by the user. The 

following matrix equations, which depend on this steering option, 

are used to compute the [IB] matrix. 

[IB] « [LB] [IL] (body rates or inertial Euler angles) \ 

[IB] - [GB][IG] (relative Euler angles) > (III-l) 

[IB] » [AB][GA][IG] (aerodynamic angles) / 

The basic relationships between the coordinate systems de- 
fined by these equations are illustrated in figure III-9* The in- 
verse transformation can generally be computed by merely trans- 
posing the matrix elements because of the orthonormality of 
these matrices. 



Figure III-9.- Matrix Transformations 
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) 

i 


» 


A summary of these matrices is given below. The symbols c 
and s denote sin and co s ,_..r.e 8f>ftC t iv el y . 


tIL] 



f 11*1 . inertial to launch.- 

The lIL] matrix depends 

on $ L . 

V 

and A^, and is given by 




C *L C0 L 

“Vi 


* 

•♦l c0 l 8A zl ‘ cA ZL 8 °L 

^ZL C °L * VA 

~ sA zl c *l 


- 8A ZL 80 L " ^Zl^L^L 

b * 21 c \ - CA ZL‘\‘\ 

cA zl c *l 


fLBl. launch to body.- The [LB] matrix is computed indi- 
rectly from the body rates by integrating the quaternion equa- 
tions, or directly from inertial Euler angles. When the body 
rate option is used, the- quaternion rate equation 


'eo 



~ e i 

e 2 

•3 

«l 


1 



-«3 

*2 


2 

e 0 

" e l 

ei 




e °. 

e i 

-e 


is integrated to compute the tLB] matrix, which is then given 

ty 


[LB] 


e 0 e l 2 
2(eie2 - 0003) 
2(eie3 + eoe2> 


2(ei«2 + e 0 e 3^ 

_2 x p 2 _ e 2 

e 0 e i *2 ®3 

2(e2V3 - eoei) 


2(eje3 - eoe2> 

2(e 0 ej + e 2 e 3 ) 
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(III-2) 


an^3> 


(III-4) 

2 

3 


• «~ 


a 


When the lnerelel Euler angle option Is used, ILBJ ' Is computed 
directly as 


ILBJ'- 


aJijCOj 




C^jSdj 


c4» I sij) I ce I + a^jSBj 


a^ jS*J< jC0 j - c^j80j 


c^spjSOj - S^jCOj 


[IQ], Inertial to geographic .- The [IG] matrix depends on 
the geocentric latitude and the inertial longitude, and Is given 
by 


[ICJ - 


- 8 * c c6 X “•♦ c #e i c * c 

-sQj c0j 0 

-C* c Cj0 -c^gSej “«* c 


fGBl . geographic to body .- The [GB] matrix depends on the 
relative Euler angles, and is given by 


[GB] - 


C V*R 


c 9 R**R 


-•9. 


•*R ce R 


Ml 


•V S R C *R ‘ C V*R *VV*R + C V*R 
C V 6 R C *R + *V*R 

hie to atmospheric relative velocity system 


«VV*R ' **R**R 


C V*R 


(ARVS).-The [GA] matrix depends on the atmospheric relative 
flight aalmuth and fllghtpath angles, and Is given by 


[GA] 


c v x a 

C V X A 

■• X A 

cX A 

*V X A 

*V X A 


- Y A 

0 

«Vj 


(IIM] 


OU-61 


an 


(in-8) 


ZZZ-9 


i 


fABl. ARVS to body .- The [ AB 1 matrix depends on the aerody- 
namic angles i and la given by 


[AB] 


cac0 -c*»a0co + aoso 

80 cdcu 

saeii -sesBco - case 


-c«s0s<i - sad 

C08O 

-SU808O + COCO 
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Other transformations, which are not related to the calcula- 
tion of the [IE ] matrix, are presented below. 


riPl. Inertial to planet relative .- The (IPj matrix trans 
forms between the Earth-centered inertial frame and the Earth- 
centered rotating frame. This matrix depends on the rotation 
rate of the planet and the total elapsed time of flight, and is 

given by 


[IP] 


eft t 

8ft t 

p 

P 

-si! t 

eft t 


L 


o 

0 


1J. 
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fRBl. body reference to body .- The [RB] matrix transforms 
data in the body reference system to the body frame. This matrix 
has a constant value and is given by 


-1 


[RB] 


0 


0 

1 


0 

0 


0 0 


-1 I • 
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IV. PLANET MODEL 


The planet model la composed of throe types of data and equa- 
tions. These are: (1) oblate planet geometry and constants, (2) 

an atmosphere model that computes atmospheric pressure, density, 
temperature, and speed of sound, and (3) a gravitational model 
that computes the gravitational accelerations. The user selects 
the appropriate models and inputs the corresponding data. The 
input data and the equations used in these models are described 
below. 


Oblate Spheroid 

The 1960 Fisher Earth model is preloaded into the program. 
This model is defined by the equatorial radius Rg, the polar 

radius R^, the rotation rate rt p , the gravitational constant 

p, and the second, third, and fourth gravitational harmonics, 

J 2 , J 3 , and J 4 , respectively. The stored values for these 

constants are: 

R e - 2.0925741 x 10 7 ft, 

Rp - 2.0855590 x 10 7 ft, 

11 - 7.29211 x 1Q “ 5 rad/s, 

P 

p - 1.4076539 x 10 16 ft 3 /s 2 , 

J 2 - 1.0823 x 10- 3 , 

J 3 » 0, 

J 4 - 0. 

The constants J 3 and J 4 are preloaded as zero, but can be ini' 
tialized by input. For example , if the Smithsonian Earth model 
la desired, then these constants would be input as 

J 2 - 1.082639 * 10 " 3 , 

J 3 - -2.565 x 10“ 6 , 

J 4 ■ -1.608 x 10" 6 , 
p - 1.407645794 x 10 16 ft 3 /* 2 , 
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Op - 7.29211515 x IQ"'' rad/s, 

Rg - 2.092566273 x IQ 7 ft, 

Rp - 2.08550242 x IQ 7 ft. 

The geometry of this spheroid Is illustrated in figure 1V-1 
The pertinent equations related to this model are 

♦« • “ ln " ('i/'tl j 

* a ■ ( k t * n •* c )' k ■ (h/*rY ( 

R g • Rg (l + (k - l)ain ? 4' c )" a 1 

h - r x - R 8 . ) 

where 4 is the geocentric latitude, <f is the geodetic latl 
c g 

tude, 6^ is the Inertial longitude, 0 is the relative longi- 
tude with respect to the planet, r^ is the distance from the 
center of the planet to the vehicle, R is the distance from 

the- center of the planet to the planet surface, and h is the 
distance from the planet surface to the vehicle. 


Worth polo 


S< *S 


Vehicle 


♦e / *1 


South pole 


Figure IV- 1 Oblate Planet 


Gravitational Model 

The gravitational model includes optionally second, third, 
and fourth harmonic terms. The potential function for this model 

Is 



The gravitational accelerationa calculated from thia potential 
function are: 

g - - BL 
U XI 

- -y ~ P (z, r) 
r J 

c - -UL 

nri 3yj 

, » -y P (z» r) 

• r J 

„ 3U 

G Z1 “ - 3^ 

. _ £_ £|i + JR* (3 - 5Z‘)) 2 + B — • (*2* - 7z 2 Z 2 - 

+ BR 4 - 10Z 2 + W 4 ) z 
where x ■ Xj, y ■ y^c * " r “ r i» * n ^ 

* - Vi 

Z - Xj/tj 

J - j J 2 

h-|j 3 

D ■ - "g" J4 

fl + JR 2 (1 - 5Z 2 ) + H (3 - 7 Z 2 )* 


P (2, r) 


+ DR 4 (9Z 4 - 6 Z 2 + 


u»|u> 



Atmosphere Models 

POST has the optional capability of three atmospheric models— 
the general table lookup, the 1962 U.S. standard atmosphere, and 
the 1963 Patrick AFB atmosphere using polynominals . The general 
table lookup model gives the user the flexibility of inputlng his 
own atmospheric model if none of the preloaded models is adequate. 
This is particularly useful In performing trajectory analysis for 
planets other than Earth. The parameters required to define the 
atmospheric effects are the atmospheric pressure p, atmospheric 
density p, speed of sound C , and atmospheric temperature 

T. These parameters are functions of the oblate altitude h. 

Table lookup atmosphere model . - The table lookup atmosphere 
model can be defined entirely by using tables that show pressure, 
temperature, speed of sound, and density as functions of altitude. 
The speed of sound and density tables can be omitted if desired} 
in .this case, the speed of sound and density are comput.ed_as . 


c 8 

P - *2 x* 


(IV-5) 


where 


Kl 


vR* . 




Kj ‘ R? 


Y “ ratio of specific heats 
Mq - molecular weight 


R* » universal gas constant. 

1962 P.8, standard atmosphere model .- The 1962 U.S. stand- 
ard atmosphere model is given as a function of geopotential alti- 
tude | H y which is computed as 


H 


8 


*A h 
a A + h* 


(IV-6) 
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The molecular scale temperature, T^, Is defined by a series 
of linear segments j as a function of geopotential altitude 

(“*)• 

The corner points connecting the straight-line segments are 
referred to as base altitudes base temperatures 

etc. From a table of base altitudes* base temperatures, and 
dT^dtt j (the slope within the linear segments), the tempera- 
ture at any desired altitude can be calculated from the following 
equation: 


[ M * \ * \ (", * "»)• 


Values of P_, T , and U, 
table 1V-1. B “b 


versus H are presented in 
o 


The atmospheric pressure* is determined as follows: 


P - P, 


vi r /«o m o\ 
W L" 1 *? 


P - P B exp - 


80^) P - «,) 

R* T„ 


for segments with 


for segments with 


i 0 , and 


where P is the base pressure corresponding to the given base 
B 

altitude H . These base pressures can be calculated once the 

g 

sea-level pressure, Pq , and the temperature profile have been 
specified. 



r 


i 


Having calculated the temperature and preaeure, the density, 
p, speed of sound-, C fl , and atmospheric viscosity, are 


determined as follows: 



i 

(IV-9) 
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Where Is the acceleration of gravity at sea level, Mq Is 

the molecular weight of air at sea level, R* Is the gas con** 
stant, y 1* the ratio of specific heats, and $ and S are 
Sutherland's constants. 


M q - 28.9644 

R. - 8.31432 « 103 ( . K) ( { - 8 -- oT) - 
Y * 1.40 

6 - l -«8 « 10-8 _ k 

S - 110. 4*K - 198.72°R 
gjj ■ 9.80665 m/sec 2 - 32.174 ft/sec 2 . 
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In the 1962 U.S. standard atmosphere, the molecular weight 
varies with altitude above approximately 90 km; in POST the molec- 
ular weight is assumed constant, resulting In a slight discrepancy 
above 90 km. In the 1962 U.S. standard atmosphere, geometric alti- 
tude Is transformed to geopotential altitude, which Is used through- 
out. Thus, above 90 km, a constant slope of molecular scale tem- 
perature versus feopotential altitude is used instead of the con- 
stant slope of temperature versus geometric altitude. 


i 
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TABLE IV-1.- 1962 U. S. STANDARD ATMOSPHERE PROFILE 


** 

* 

rr 

T> B» P 8 f 

V 0,5 

hi * °R/ft 
B 

0.0 

0.21162166 + 4 

518.67 

-0.35661600 - 2 

36 089.239 

0.47268050 + 3 

389.97 

0.0 

65 616.797 

0-.11434543 f 3 

389.97 

0.54863995 - 3 

104 986.87 

0.18128948 + 2 

411.57 

0.15361920 - 2 

154 199.48 

0.23163263 + 1 

487 . 1 7 

0.0 

170 603.68 

0.12322603 + 1 

487.17 

-0.10972801 - 2 

200 131.23 

0.38032532 + 0 

454.77 

-0.21945600 - 2 

259 186.35 

0.21673064 - 1 

325.17 

0.0 • 

291 151.57 

0.34333824 - 2 

325.17 

0.16953850 - 2 

323- 002.74 

0.62hi4785 - 3 

379.17 

0.28345707 - 2 

354 753.59 

0.15361733 - 3 

469.17 

0.56867005 - 2 

396 406.39 

0.52676024 4 

649.17 

0.11443751 - 1 

480 781.04 

0.10566108 - 4 

1 729.17 

0.86358208 - 2 

512 046.16 

0.77263469 - 5 

1 999.17 

0.57749093 - 2 

543 215.48 

0.58405376 -5 

2 179.17 

0.40610461 - 2 

605 268.45 

0.35746030 - 5 

2 431.17 

0.29274135 - 2 

728 243.91 

0.14559124 - 5 

2 791.17 

0.23812804 - 2 

939 894.74 

0.39418091 - 6 

3 295.17 

0.20152600 - 2 

1 234 645.7 

0.84380249 - 7 

3 889.17 

0.16354849 - 2 

1 520 799.4 

0.22945543 - 7 

4 357.17 

0.11010085 - 2 

l 798 726.4 

0.72259271 - 8 

4 663.17 

0.73319725 - 3 

2 068 776,5 

0.24958752 - 8 

4 861.17 

0.0 





1963 Patrlek AFB atmosphere usin g polynomials.- In this 
niodil, pressure end temp erst ure ere calculsted as functions of 
geometric altitude (h) . These persmeters ere calculated in met- 
ric units and converted to English units if required. 

Pressure: 

1) Altitude region » 0 to 28 000 meters: 

P » Pj exp (A + Ai h + Aj h 2 + A 3 h 3 + A^ h 4 + A 5 h 3 ) 
where Pi » 10.0 Newtons/cm 2 ; 

2) Altitude region » 28 000 to 83 004 meters: 

p » gg x 10 " 4 exp (A + Ai h + A 2 h 2 4* A 3 h 3 + A 4 h 4 + A 5 h 5 ) ; 

3) Altitude region » 83 004 to 90 000 meters: 


-1.373301523 x 10 12 h - h 


P * P B exp ll B (6344860 + h) (6344860 


‘ h B. _y 

+ w 


KiV-ll) 


4) Altitude region » 90 000 to 700 000 meters: 

. 1.373301523 x 10 * 2 \ 

L n (P) " L n ( P BI + L m (6344860 + h) (6344860 + h B ) J 




/h - h. 
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Temperature 


1) Altitude region « 0 to 10 832.1 meters: 

T ■ T* • A + Aj h + A;j h ? + A 3 h 3 + At, h 4 + A r , h 5 ; 

2) Altitude region = 10 832.1 to 83 004 meters:* 

T « A + Ai h + A> h' + A 3 h ’ + A 4 h 4 + A 5 h 5 ; 

3) Altitude region * 83 004 to 90 000 meters: 

1 - T B + he (*> - M- y 

However, In this region * 0, and thus 
T - - 180.65°K; • 

4) Altitude region • 90 000 to 700 000 meters: 

Density: 

1) Altitude region * 0 to 28 000 meters: 

o » Pi exp (A + Aj h + A; h** + A 3 h^ + At, h 14 + A 13 h^*) ; 

’ 

2) Altitude region = 28 000 to 700 000 meters: * 

p - (34.83676) |. 


♦Virtual temperature is the same as kinetic temperature 
above the 10 832.1>meter altitude. 


(IV-12) 
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TAing XV-2 DERIVED COEFFICIENTS FOR THE 1963 PATRICK AFB ATMOSPHERE MODEL 
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Pressure and density ratios: 

Altitude region - 0 to 700 000 meters: 



Velocity of sound: 

V_ - (20.046707) (T)**. 
s 


(IV-14) 


(IV-15) 


The atmosphere model-derived coefficients are presented in 
'.able iV-2. The molecular temperature gradient is documented in 
table IV- 3 for geometric altitudes from 90 to 700 km. 


Winds 


The atmospheric wind velocity components are input in tables 
using either meteorological or vector uotation. If these tables, 
which are normally functions of oblate altitude, are not input, 
then the atmosphere is assumed to rotate uniformly with the 
planet . 

The wind velocity components can be input directly in the 
geographic frame by defining u^, v w » and v w , or by defining 

the wind speed ,,2 i mu ^h an< ^ t ^ ie w * n< ^ 

azimuth bias |A ZWB j . The resulting wind velocity components in 

the 6-frame are: 


V 

•^WG 


V M (h) COB (h) + Agygj 

V W (h) 8in (*ZW + ^zwb) 


L w w (h) 


(IV-16) 


IV-14 


It Is clear from the above equation that In order to Input vector 
wind data A ZMB must be input as zero, whereas for meteorologlc 

data the preloaded value of 160 ° should be used. 


The wind velocity In the GCI frame Is then given by 



liu]-' 



Thus, the atmospheric relative velocity vector in the ECI frame 
Is 



and Its magnitude is given by 


*1 



(IV-17) 


(IV-18) 


V. VEHICLE MODEL 



The various physical properties of Che vehicle are modeled 
by the- user when he selects the pertinent options from the set of 
vehicle simulation modules. The equations used in these modules 
are presented below. 


Hass Properties Model 

The gross weight of the vehicle at the beginning of each 
phase Is given by 

W - Vf + w , 

G stg pld* 

where W gt g Is gross weight without payload and is the pc 

load weight. For phases other than the first, the gross weight 
can optionally be computed as 

u + ■ u" - w _ u 

G *G w jett W PR* 

where W* Is the gross weight on the positive side of the cur- 

rent event, U Q is the gross. weight on the negative side of the 

current event, W^ fttt Is the jettison weight, and W pR Is the 

weight of propellant remaining. These options are obtained au- 
tomatically, based on user Input. 

The propellant remaining Is given by 

U ■ U • U 
PR w P t PC’ 


where W p is the initial weight of propellant and W pc Is the 
amount of propellant consumed. This latter term Is given by 


“pc ■ f* 


dt + W, 


where W is the total rata of change of the vehicle's weight. 


V-l 


The composite Inertia matrix Is Input with respect to the 
body axis system which 1« located at the instantaneous composite 
center of gravity of the vehicle.* In 6D POST, the moments and 


products of 

Inertia are < 

ilnf tned 

ha tho Intograla 


I • (y‘* + 

XX J 

z : ‘ dv 


1 xy " 

J xy dv 


'yy ' I*' * 

z- dv 


! xz " 

$ xz dv 

(V-5) 

* 

y*‘ dv 


l y* " 

$ yz dv. 

(V-5) 

The inert la matrix 

is then 

given by 





"1 

XX 

" l xy 

-xl 




II] - 

-t 

xy 

‘yy 

_1 yz 


(V-O 



-1 

xz 

_I yz 

1 

I 

zz 




and la generally input as a function of the wei ght of the vehicle. 


The composite center of gravity la refetenced with respect 
to the vehicle reference Itarne, and the compoients |x , y » * ' 

are generally Input as a function of vehicle ueight.' ® 


Propulsion Calculations 


6!) POST can simulate both rocket and jet engines. Th& pro- 
gram can simulate up to 15 engines in either mode. 

Rocket engines .- There are two input options for engine data 
in the rocket mode. In the first option, tables for vacuum thrust 
and maximum weight flowrate aro Input for each engine. In the 
second option, tables for vacuum thrust are input, along with the 
vacuum specific impulse for each engine. The vacuum specific in*’ 
pulse is then used to calculate the mass flowrate. 


The rocket thrust par engine is given by 


T 


R 


i 


>' T 

vac 


- A_ p(h) , 
b i 


(V-7) 


*In 6D POST the center of gravity is assumed to coincide 
with the center of mass. 


V-2 


where n Is the throttle setting, T vac Is the vacuum thrust 

of the 1 th engine, Is the nozzle exit area, and p(h) Is the 

atmospheric pressure. Summing over all engines yields- the total' 
rocket thrust 


T 


R 



<V-8) 


where N is the number of thrusting engines, and N ^ 15. 

6Ug eng 

The weight flowrate in the rocket mode is given by 


eng 


"E (“T 4 


1-1 

*eng 


-n 


(V-9) 


E ( T »»c/ I »p vac ) 1 


Jet engines . - In the jet engine mode the net jet thrust per 
engine is given by 



f(H, n), 


(V-1-0) 


where 

6 - P(h) /p gL 


V-3 


and t=- (M) Is a monovariant table. The total Jet thrust is then 

6 


given by 


N 


eng 


Tj ■ '■ Z/ <h) /( p «.) ( t Vm 

1 * 


i-1 


The weight flowrate in the Jet engine mode is 


N 


eng 


‘■-SJsFwiniH 


6 . 


The thrust equation for each engine is given by 

cos <Se cos 




l 


sin 


- cos sin 


where the pitch and yaw engine deflection angles 6e and 5e^ 
are defined in Figure V-l . 




Aerodynamic Calculatlcr.j 


The aerodynamic force coefficients can be expressed in terms 
of the lift, drag, and side-force coefficients C L , Cg, and 

Cy (fig. V-2) , where and are directed normal to, and 

along the velocity projection in, the x b" 2 b plane. Note that 

Cy produces a side-force, F^ , acting in the direction of y^. 

Lift and drag force coefficients are transformed to axial 
and normal force coefficients as follows: 



where a is the angle-of-attack. 



Figure V-2.- Aerodynamic Angles 

Xhe aerodynamic coefficients can also be expressed in terms 
of the axial force, normal force, and side force, C^, C^, and 

Cy, respectively. Here C A and produce forces that act in 

the -Xg and -* B directions, and Cy produces a force acting 

along y.. 






Each aerodynamic coefficient is computed by Interpolating 
the values in the table. In general, eight tables are allocated 
to each coefficient. These tables can be- monovariant, bivariant, 
or trivariant, and seven tables per coefficient can have arbi- 
trary hollerith mnemonic multipliers. This generality enables 
all standard forms of aerodynamic data to be directly input into 
the program. 

The aerodynamic force coefficients are obtained by summing 
the individual contributions as follows: 


0. B C A + C a (*,M) + C 6a + C 
A A 0 A . A 6a A «e 


3 

<5e + C A C. 6f , (V-15) 

or 1 4^ A 6f, 1 


i-1 i 
3 


(V-16) 


C K - % + S la,M) + V 48 + C N. 68 + V 48 + 2X f “i* 

l) i** Or i 

or optionally 

3 

S - C D„ + + V„ 48 + 48 8 V 48 + & Sfi “f <V ' 


17) 


C L * C L. 


+ C T (a ,H) + C. 6a + C <Se + C 6r + / (V-18) 

0 L Ha He Hr Hf ± 1 


and 


CL, - C v + C v (e,M) + C 6a + C «e + C„ 

Y Y 0 Y Ha He 


oa 


6r + > C v 6f.. (V-19) 
48 4£ i 


The aerodynamic moment coefficients are given by: 

6a '6o "dr 


C l m C i + c t< p » M) + <V 4a + c i,_ 6e + C i, 6r 


*£■ 


(V-20) 


61 + C. r* + C 

w 4£ i 1 8 ' 


V P 


Vr 6 


C • C + C (o,M) + C 6a + -C 6e + C 6r 
m ® 0 ® v ®6e ®6r 


+ 2X, st i- + -v. 


s— » “6f * q ’ 

1-1 r i q 

C - C + C (8,M) + C 6a + C 6e + C 6r 
» n o n ' n 6a n 6e n 6r 


+ V C 6f + C p’ + C **» 

AJ n» f i n . n r t 

i-1 6f i p r 

where i - 1,2,3, are arbitrary user defined deflection 

angles; and 

*’ • ,,d R/ 2V * ’ “.V”* 

,• - ,d p/ 2v A - y P /2v A 

*’ - r V 2V A - uA / 2V A- 

The Mach number and dynamic pressure are given by: 


,3 

c s 


’ • 2 * V A' 


where p is the atmospheric density, V A is the velocity of the 
vehicle with respect to the atmosphere, and C g is the speed of 

sound* These atmospheric parameters are determined from the atmo- 
spheric models as a function of the altitude h above the oblate 
spheroid; l.e.. 


The angle of attack In pitch (a) and the angle of sideslip 
(8) required to determine the aerodynamic coefficients are 
calculated as follows: 


a ■ tan' 


8 • tan' 


-i r-M 

Lcos al 

1 fsln e l 
[_cos ej* 


W, 


sin a 


cos » 


B 


sin 0 » — 
V A 


+ w b 

% 

/u B + w B 
B 


cos e 


H +> w2 


B 


V. 


The total angle of attack is 
<* T - cos 




The aerodynamic forces in the body frame are 


*AB * * S 


- C A 

°Y 


(V-25) 


(V-26) 


(V- 27) 


where q is the dynamic pressure and S is the refer once area* 


Aerohcatlng Calculations. 

6D POST provides for a wide variety of aeroheatlng calculations. 
Some of these options are specific in nature and apply only to 
particular vehicles, whereas others are quite general. The gen- 
eral heat rate option is based on trivariant table interpolation 
and provides complete flexibility with regards to vehicle shape 
and heat-transfer methodology. The various heat rate equations 
are described below. 


V-8 


Heat rate equation*. - 


1) 


Chapman’s equations . In this calculation the heat rate 
is given by 



(V-28) 


where is the nose radius » p is the atmospheric density, 

and V r is the reference circular orbital velocity. 

w 

2) General table lookup . This heat rate is given by 

Q - Q fc (xi, x 2 , x 3 ), (V-29) 


where x\ ¥ x 2 , and x 3 can be any internally computed variables. 

For example, the values that would normally be selected are *1 * 
x 2 » h, and x 3 • Vj^ 

3) Modified Chapman’s equation . Here the heat rate is given 

Q - Q t (x-i, x 2 , x 3 ) Q c , 

where- Q fc is an arbitary table and is the standard Chapman s 

equation. 

4) Turbulent-flow heat rate . The turbulent-flow heat rate 
is given by 


T / \ 0,8 / V a\ 318 

Q - Q t («!. *2. «3> U*>° (4) j- 0 

5) Maxiwa «*«*• *11— heating . The equations for this method 
are ‘given, below In eequence. 

a) Altitude-velocity correction: 
lh . 10* f 1.06112 - 6.16586 V^O'J 

+ 51.12MO (V^O*)* - 20.66256 (» A /l0 4 ) 5 ] I ( 

+ 22.52588 (v^0"j 2 - *8.280Bo( I 


h t.f * " + 4fc - 


J 

V-9 



b) Maximum centerline heat rate at reference conditions: 

— If h ref > 103 600 m: 

q ref - 10 2 [277.93332 + 134.55760 h re£ /10 5 - 807.75941 (h ra ^10 5 ) 2 

+ 2.90536 (h ref /l0 5 ) 3 + 722.36896 (l» ref /10 5 ) 4 - 311.40176 

(h f /10?) 5 1 ; 

— If h ref < 103 600 m; ' W ' J 

q f - 10 4 [7115.39692 - 34 881.13588 h ^ 10 5 + 69 844.23141 

( h refA-° 5 ) 2 

— 71 534.98453 (^ ref /10 5 ) 3 + 37 506.13054 (h ref /X0 5 ) 4 - 

- 8048.55112 |h re ^10 5 ) 5 J. 

c) Angle of attack correction: 

q., a „ Jh-*- • l £ n (x)] 2 , 

laax » <xf max » a* 50 

where 

x - 10 2 [0.01136 + 0.01343 a/10 2 + 1.42672 (o/10 2 ) 4 - 0.75623 

(a/10 2 ) 5 ] + 0.30535 (a/10 2 ) 2 - 1.06269 (a/10 2 ) 3 . 

d) Maximum centerline heat rate: 


>(V-33) 


(V-34) 


'max 


(^max,a) / ( < *mex,a«50 # ) (^ref )• 


(V-35) 


In addition to the heat rate calculations, the program also 
provides the capability to calculate other aeroheatlng Indicators 
that can be used for trajectory shaping purposes. 

Aerodynamic heating Indicators .- the heating rate for rero 
total angle of attack is 

Q - q V A . (V-36) 


V-10 


The aerodynamic heat lng_lndlca tor for sero total angle of attack 
la 


Q ■ J* Q dt. 

0 

The heating Indicator for non-zero angles of attack is 
given by 


where 


and 


b 

Q' - f t (o', M) Q dt, 


/ 7 \3/7 

% M) « f 1 + y M 2 sin 2 o' j K, 

-M'-H* f’|: 


K 




% 


% 


a 

Q» 


- Q' 


for o < 0° 


- o | 

> for o > 0° 

- Q" | 

- e ) 

> for 6 < 0* 

- r ) 

-e | 

- r ) 


for 6 > 0*. 


(V-37) 


(V-38) 


> 


(V-39) 


The-heating Indicator for laminar flow is calculated aa 


t 

• / 17 

0 

where 


The heating Indicator for turbulent flow is calculated as 

r /n \ 0 - 8 / V 4 \ 3 * 18 

<>turb -/ 1500 % ( ITooo) « 

Ten-Panel Vehicle Heating Model .- Special aeroheating calcu- 
lations are available for a ten-panel vehicle model. The heating 
ratios are referenced to the heat rate calculation. The total 
heat for each panel Is given by 



where Q is the total heat and H. la the heat ratio for panel 

R 1 

1. The weight for each panel Is the product of the weight per 
unit area and the area of the panel. The total weight is the sum 
of the individual weights for each panel: 

10 

«p ■ £ v A t 

l-i 1 

where is the weight per unit area and la the area 

of the l** 1 panel. 



(V-40) 


(V-41) 


(V-A2) 


(V-43) 


(V-4A) 


Sensor Module 


The sensor module computes information that describes the 
behavior of the sensing elements of the vehicle's navigation sys- 
tem.- Thus, the primary functional responsibility of the module 
is that of simulating hardware characteristics of For 

example, the behavior jo£ an inertial measurement unit (IMU) can 
be described by a mathematical model of the platform and the ac- 
celerometers. Frequently this module is used for error analysis 
purposes . 

Sensor models called by this module are necessarily vehicle 
and subsystem dependent. As a consequence, the sensor model must 
te designed and implemented for each particular application. 

There are many applications of the program that do not re- 
quire a specific simulation of the sensors. Therefore, for con** 
venience, a "perfect” sensor model is coded into this routine. 
This "perfect" sensor model sets -the sensed program variables 
equal to their actual values as calculated in the simulation 
models . 


Navigation Module- 

The function of the navigation model is to estimate the state 
(position, velocity, etc) of the vehicle based on the sensor out- 
puts. Clearly, this module is also vehicle and subsystem depen- 
dent and must be designed' and implemented for each specific ap- 
plication. This version of the program contains no navigation 
models. As a consequence, the estimated state is set equal to 
the actual state. This is equivalent to simulation of perfect 
navigation. 


Guidance Module 

The guidance module takes the output of the navigation model 
and computes a guidance command. Typically, the guidance corn- 
ed represents a desired change in the current attitude of the 
vehicle. This command is computed on the basis of meeting some 
specified trajectory condition, such as, inject conditions or 
landing conditions. The autopilot is designed to remove the er- 
rors between the commanded values of the guidance variables and 
Slir eetuel <ot ..n..d> valu.a . Thl. t. ac,o.pU.h.d by deflect- 
ing engines, control surfaces, and/or firing RCS jets. 


V-13 


I 


I 


The current version of 60 POST contains three preloaded 
guidance options# (1) an open-loop profile steering} (2) a closed 
loop v-h profile ascent algorithm; and (3) the constant drag Space 
Shuttle reentry scheme (ref# H-l) . If these methods are Inade- 
quate 9 the user may implement his own guidance algorithm into this 
module# 


Autopilot Module 

The function of the autopilot module is the generation of a 
command, which, when implemented through the deflection equations 
contained in the controls module, causes the vehicle to respond 
as prescribed by the guidance module# This functional responsi- 
bility is depicted in Figure V-3# 


Guidance module 

Mndslg ♦ 

9 « - 6 - 0 
-t Ti 

Autopilot module 
Models: 

66 

Control* model 
i rn 6 ♦ (M) 49 

1 

1) V v# h profile 

2) Shuttle reentry 

— 

r 

0 

—a 

1) Shuttle ascent 

2) Shuttle reentry 



Actual or tensed- 
vahlcla state, e.g*, 

acceleration* attitude, 
attitude rate 



.from Simulator 


Nomenclature : 

e - Guidance command 
— c 

0 - Actual or sensed values of guidance variables 

—a 

e - Generalize error signal 

69 - Pitch, yaw, and roll autopilot command 

6 •• Deflection (engine or aerodynamic surfaces) angles 


V-14 


Figure V-3.- Functional Flow 







The autopilot module calculator! only autopilot commands 
based on the input guidance commands. And does not calculate en- 
gine or control surface deflections. The engine and control sur- 
face deflections are computed in the controls module as a linear 
function of the autopilot commends. The autopilot commands 66 , 

«<; , 641 represent changes in vehicle attitude. The mixing equa- 

tions determine the engine and control surface deflections that 
create the control forces and moments. 

Currently, there are two Space Shuttle autopilot models avail- 
able in 6 D POST. One autopilot is for ascent and the other for 
reentry. The ascent autopilot is somewhat standard and could be 
used on most ascent problems with little or no modification. The 
basic inputs to this model are: attitude commands from the guid- 

ance, inertial attitude angles, body rotational rates, transla- 
tional accelerations, and preloaded engine deflection commands. 

The outputs are pitch, yaw, and roll autopilot commands, which 
are sent to the controls module to determine the engine deflec- 
tion angles. The reentry autopilot is Space Shuttle oriented and 
is probably not applicable to other vehicle configurations. This 
model is intended to provide attitude control for Space Shuttle 
beginning at approximately 400,000-ft altitude and ending In the 
high subsonic flight regime. The control logic makes use of both 
aerodynamic control surface torques and reaction control jets. A 
complete description of this model is presented in ref. H-2. 


Controls Model 

The controls model converts pitch, yaw, and roll autopilot 
commands into aerodynamic control surface deflection angles and/ 
or engine global angles. The conversion of the autopiloc com- 
mands Into deflection angles Is Implemented through the matrix 
mixing logic given by the equation 

6 - «jo + [M] 66 ., (V-45) 

where £ denotes e general deflection angle with a null posi- 
tion of 6 jq, [M] the mixing gains, and 6 £ the autopilot com- 
mands. The gains contained in the mixing matrix, [M], and the 
null deflections, £ 0 * &re specified by user input. 
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The standard aerodynamic surface deflection mixing equations 
used In the program are 

<$a - 4 q -t RR 6j| 4 + KP fij| 4 + K* 6a 4. (V-46) 

4 - 6e o + KR 6o ii + KP 6<| 60 + KY^ 4 (V-47) 

6R • 6R + KR. 4 + KP . 60 + KY. 6*, (V-48) 

o dr 6* 6r 

and the standard engine deflection angle mixing equations are 
similarly 


6a - 6e + KP. 60 + KR, 6$, (V-49) 

Pi P l0 6p 6p 

68y i " 6Cy i + KYfiy + l * 6 y 6 *’ (v ’ 50) 

where 60 , 64 *, and 4 are the pitch, yaw, and roll autopilot 

commands . 


Airframe Modal 

The airframe model computes the total thrust and moments 
acting on the vehicle. The forces and moments are computed from 
the engine and aerodynamic deflection angles and the RCS thrust 
and moments. 

The nongravi.tational force acting on the vehicle is computed 
in the body frame as 


h 


-TB + -AB + -RCS’ 


(V-51) 


where is the total force due to the engines, is the 

total force due to aerodynamic effects, and F^g is the total 

force resulting from the reaction-control system. Similarly, 
the total moment acting on the vehicle is computed as 

4 • % + 4. + Bros' 

The thrust vector components for both rocket and Jet engines 
are determined from the thrust magnitude T_ or T. and the 

R i J i 

engine gimbal angles 6e and 6e . The total thrust force 

p i y i 

in the body frame is given by: 
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* E v ■ 

i-i 1 


where the individual engine components, , are given by 


F 

•^TB, 


cos 6e cos 6e 


sin 6e 


-cos 6c sin 6e 


For roll nozzles, the thrust vector is given by 


‘ \ 


cos 6e cos 6e 


sin 6e 


-cos 6e sin 6e 


where for a roll nozzle the deflection matrix [K] is given by 

1 0 0 * 


K 


0 c<|> i sd> t 

0 -84^ c^ 


where 


or 


- Input value 


- tan- 


1 (W , * p i) u “ ot in,,ut ' 


(V-53) 


(V-54) 


(V-55) 


(V-56) 
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The thrust moments are obtained by summing the thrust moments 
for each engine as follows: 

N 


eng 

" -£ ht, x 

i = l 


(y-57) 


where 


ar, 


BT, 


gpi 

- v 1 
cgj 

r 

- y J 

8Pi 

c 8j 


- 2 

! gPi 

eg 


The aerodynamic forces and moments are given by: 


and 



"p ” 

AXB 


— * 

- C A 

F * 
-AB 

F 

AYB 

* qS 

°Y 


F 



! 

AZB 


N 


-AB 


qS 


d R 


d u c 

P m 


d Y C n 


AB X A -AB 


(V-58) 


(V-59) 


(V-60) 


where 


A -AB 


~( X ref " X cg) 

) 


( y ref ^cg 
~ ( z ref Z cg) 

The aerodynamic reference point ( x yef’ ^ref* Z ref) 
calculated from tabular Input. 


is 


(V-61) 


The RCS forces and moments are computed in the autopilot 
model and are merely added to obtain the resultant force and mo- 
ment vectors. 
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VI. TRAJECTORY SIMULATION 

The following sections present the equations used in the 
trajectory simulation subroutines. Theaa equations summarise the 
principal computations performed by the program, and motivate many 
of the program input procedures. 

Events/Phases 

Simulation data are input according to phase, where the 
phases are defined by a user-specified sequence of events. The 
simulation equations are then solved sequentially by phase. 
Therefore, the user is required to input a sequence of trajectory 
segments that define the problem being simulated from beginning 
to end. These trajectory segments, or phases, are defined by 
two events — a beginning event and an ending event. An event is 
an interruption of the trajectory simulation that occurs when a 
user-specified variable reaches a user-specified value. An event 
must be created whenever the user wishes to change any input data 
for the problem or to cause any change in the method of simulating 
the problem. For example, the sequence of events for a typical 
ascent problem could result in a simulation setup similar to that 
shown in figure VL-1. 


■ Phase 2 


4 


Pea criotlon 
liftoff 

Initiate pitch rat* 1 at 20 a*c 

initiate pitch rat* 2 at 30 aec 

Initiate pitch rat* 3 at 60 sec 

Initiate anti* of attack control at 75 aac 

Jettison ttaga 1 whan propellant consueed 

Initiate pitch rata 4 20 eac after staging 

Initiate yaw rat* t 100 aac after event 7 

Orbit injection at inettial velocity of 25 568.0 fpe 


Figure VI-1.- Event Sequence Setup 
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The event numbers for a given problem must be specified as 
reel numbers by the user in monotonic Increasing order. These 
event numbers are then used by the program to determine the order 
in which the events are to occur. The program requires that each 
problem have a minimum of two events—an Initial event and a final 
event. Since a phase is initiated by the corresponding event, the 
event criterion for a given event specifies the conditions at the 
beginning of the corresponding phase. A problem is terminated 
by specifying the last event that is to occur. The problem can 
also be terminated in a psuedo-abort mode by specifying the maxi- 
mum trajectory time, maximum altitude, or minimum altitude. 

Although event numbers must be monotonic increasing, they 
need not be consecutive. This allows the user to easily add or 
delete events from an input deck. 

Three types of events have been defined to provide flexibility 
-in setting up a given problem: 

1) Primary events - These describe the main sequential 
events of the trajectory being simulated. These events 
must occur, and must occur in ascending order according 
to the event number. Most problems will usually be 
simulated by a series of primary events; 

2) Secondary events - These are events that may or may not 
occur during the specified trajectory segment. Secondary 
events must occur in ascending order during the interval 
bounded by the primary events. The occurrence of a 
primary event will nullify the secondary events associated 
with the previous primary event if they have not already 
occurred; 

3) Roving primary events - These events can occur any time 
after the occurrence of all primary events with smaller 
event numbers. They can be used to interrupt the tra- 
jectory on the specified criterion regardless of the 
state of the trajectory or vehicle. 

The program monitors as many as ten events at a time, depend- 
ing on the types of events to determine which event is to occur 
next. This gives the user a powerful tool for simulating complex 
problems . 
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Multiple events are monitored in the following sequence: 

1) The next primary event is monitored; 

2) As many as nine primary roving events are then monitored* 
provided there are no secondary events . A roving pri- 
mary event is added to the- list of those being monitored 
as soon as the primary event immediately preceding that 
roving event has occurred; 

3) Next, as many as nine secondary events are monitored, 
provided there are no primary roving events. (Note that 
caution must be exercised when using secondary events 
because of their nature. Since as many as nine sec- 
ondary events are monitored at a time, any one of those 
nine will occur as soon as its criterion has been met. 

Because they are secondary events, the event that occurs 
will cancel all secondary events with smaller event 
numbers . ) ; 

4) Finally, a total of nine primary roving and secondary 
events are monitored. 

Since the program can only monitor nine events (in addition 
to the next primary event) , the sum of the primary roving events 
and the secondary events must be less than or equal to nine or a 
fatal error will result. 

The time-to-go model (TG0M) determines when the events 
occur during the trajectory simulation. Basically, TG0M checks 
the values of the criterion being monitored at each integration 
step. If none of the criterion values has bracketed the desired 
cutoff value, then another integration step is taken. If a 
criterion variable is bracketed with the input step size, then 
TG0M computes a new stepslze equal to the predicted time-to-go. 

The predicted time-to-go for each event is computed from the 
equation 

At* - - y 2 (t)/(y(t + At) - y(t)) (Vl-1) 

where y(t) is the difference between the actual and the desired 
value of the event criterion. If more than one event is bracketed, 
then the minimum predicted time-to-go is used as the integration 
stepslze. This process is repeated until the criterion value is 
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within the specified tolerance of the desired value. If the 
desired condition cannot be achieved in 20 iterations,.. an error 
message is printed and the program stops. Generally this situa 
tion is caused by an input. error. The fundamental features of 
the time-to-go logic are shown in figure VI*2*~ . 


Figure VI-2.- Illustration of Time-to-Go Logic 


Translational Equations 


The translational equations of motion are solved in the 
planet-centered inertial coordinate system. These equations are 


it ' t"!- 1 [ire + *ab] + %• 


where ^TB 18 the thrust acceleration in the body frame, A^g 
is the aerodynamic acceleration in the body frame, and is 
the gravitational acceleration in the ECI frame. 


“5L 



Initialization *- There are five options for initializing 
the velocity vector and two options for initializing the position 
vector.*-- -These options are described below. .... 


Inertial position components y^, z i)*" The inertiaJ — 

position components can be input directly since no transformation 
is required. 


Earth-relative position ^0^. or 6, ^ or 4 g * h or rj.- In this 

option the equations vary and the sequence of calculation varies 
according to the choice of input. However, the basic equations 
used are: 

0j ■ 0 + |t - tgj if 0 is input, 

4 - tan -1 |k 2 tan 4 ) if 4„ is input, ) 

c » g' g f 


r i • h + R . (♦«) 


if h is input. 


J 


and 


-I " r I 


cos-4 cos 0 T 

C 1 


cos 4 sin 0, 
c I 

sin 4. 


Inertial velocity components (V^, ? Y1 , V 8I ).- These variables 
can be input directly. 

Inertial local horizontal (Vj, Yj. A Z ij’~ The inertial com- 

ponents in the horizontal frame are first transformed to the 
geographic frame as 

re. Vj or- Ajj-I 



cos Yj tin A ZI 


L-sin Yj 


(VI-4) 


(VI-5) 


(VI-6) 
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and then transformed to the ECI frame by 

h * W* — IG’ 

Earth-relative local horizontal (V^, y r * A^ R J . - The. Earth- 

relative velocity components are first transformed to the geo- 
graphic frame as 



'R 


■cos y r cos * 2r ' 
cos y r sin ^ 
,-sln y r 
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and then transformed to the ECI frame by 




[IG] 


-1 


•HRG 


+ Q x r 


-I 


Atmospheric relative local horizontal ( V^, Y^» A^j.- The 

atmospheric relative velocity components are first transformed 
to the geographic frame as 



cos y. cos 
A 

a za" 


v 

WXG 

V = V 
-AG A 

cos y. sin 
A 

a za 


V WYG 


-sin y a 
A 



V V&G 


and then transformed to the ECI frame by 
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Vj - [IG]“1 v AG + ^ X £ x (VI-U) 

Orbital parameters (h^, h ft , i, 0, 9).- This option Initializes 

both position and velocity. The equations used to transform the 
orbital parameter to the ECI position and velocity are: 





Rotational Equations 


The rotational equations of nation are solved in the body- 
centered coordinate system. These equations are 

i-i IE] % 

- [ir l [M B - [i] u> B - Y [I] -B ] * 

-B " -TB + -AB + M -RCS 


(VI-14) 

(VI-15) 


where e is a four dimensional vector of quaternion parameters, 

(E] is the quaternion matrix, is the inertial angular velocity 

expressed in the body frame* M fi is the total external moment- act- 
ing in the vehicle as a result of the thrust, the RCS, and the 
aerodynamic forces, and [I] is the inertia matrix for the compos- 
ite vehicle. The (X) and [E] matrices are given by 


[ 1 ] 


L XX ■ X XY ' 1 

L xz 

-I XY I yY - 

r YZ 

~ X XZ “ X YZ 

1 a 

1— 


Fei e 2 e “I 

! 
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eo -ei 
e 0 . ei -e. 

The body rates are defined below and illustrated in Figure 
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u - Roll body rate. The angular 
x rate about the Xg-axis in 
deg/sec, 

w t - Pitch body rate. The angular 
y rate about the y -axis in deg/ 

D 

sec, 

o - Yaw body rate. The angular rate 
z about the Zg-axis in deg/sec. 


__-3" 

<[^ - 


t 


B 

Figure VI-3.- Body Rates 


Initialisation .- The rotational aquations of motion are ini- 
t ial ized by defining both the attitude angles and the attitude 
ratee. There are three optlone for Initializing the attitude! 

(1) Inertial Euler angles; (2) relative Euler angles; and (3) 
aerodynamic angles. There are three options for initializing 
the attitude ratess (1) body rates; (2) Inertial Euler angle... - 
rates and (3) relative Euler angles ratea. The attitude angles 
are used to compute initial values of the quaternions In order 
to Initialize equation VI-14. The rates are used to Initialize 
the moment equations. The equations for these options are pre- 
sented below. 


When inertial Euler angles are input the initial quater- 
nion vector Is given 

So ■ S (*!> * s (*t) * a ( e x> 

where the asterisk denotes quaternion, multiplication and where 


(♦ x ) “ cos 

(0.5 ♦,) 

+ sin (0.5 

1 


(♦j) - cos 

(0.5*,) 

+ sin (0.5 ♦jj 

k 


( 9 l) ■ “• 

(0.5 .,) 

| + sin (0.5 

j. 

J 


(VI-18) 
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When aerodynamic angles are Input, then the initial quater- 
nion vector is given by 

So ■ S (V) * s<»°> * S (* L ) * A (-« L ) * S (8 t ) * a (-♦„) * a<-w * 

® ( X A ) * • (y a ) * * *<"&) * •<“) » (VI-20) 
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where 


~ (^zl) " cos ^t) ** s * n ^ 

e(90) ■ nos (45) + sin (45) J 

e (4 L j ■ cos (0.5 4^ + sin (0.5 <J> L ) j 

e (-OjJ * cos (0.5 0jJ - sin (0.5 <’j J k 


\ 


MM 


COS 


(0.5 0j) + sin (0.5 Qj) k 


- Vc 


■ cos (0.5 4 C J - sin (0.5 4 c j J 

e(o) * cos 0.5o + (sin 0.5o) 1 

e(-B) » cos 0.5B - (sin 0.50) k- 

e(o) » cos 0.5a + (sin 0.5a) J. 

“ / 
When relative Euler angles are input, then the initial 
nlon vector is given 

«0 ' S (*21.) * sW) * e ,= L ) * . (-8 l ) * . («, ).« . (-» c ) 

* e(v) * e(e) * e ($ ) , 

where 

S (* zl ) ■ cos (°- 5 * Zl ) • 8ln (°- 5 * zl ) k 

e(90) » cos (45) + sin (45) j 

e ($ L j • cos (0.5 4 L j + sin (0.5 <t> L ) j 

e (-0jJ 9 COB (0*5 ^l) ” sin 9^) k 

e (0jj - cos (0.5 0j| + sin (0.5 0j) k 

e (4 c j - cos (0.5 4 c j - sin (0.5 4> c j j 

e(4) ■ cos 0.54 + (sin 0.54) h 

e(0) ■ cos 0.50 + (sin 0.50) J 

e(9) ■ cos 0,54 + (sin 0.54) 1 

/ 


(VI-21) 


quater- 

* e(-90) * 

(VI-22) 


(VI-23) 
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The available options for Initialising the moment equation 


1) Input ui^, Wy, directly 


ate 


2) Input the Inertial Suler angle rate® 4^, and cal- 


culate t*> , b> # u) via 
x' y 2 


w x $j co * cos ®j - sin Gj 


w y • ®i - sin ♦j 


W 1 I ^1 COtt ^1 ®! * ^1 C08 ®I 


3) Input the relative Euler angle rates 9 R , $ R and 

calculate w , w , U 2 via 
x* y * 


“* *R " 8ln ®R *R 

u - a + cos ♦_ 0_ + eln cos 6_ 4», 


R *R I * 


cos $ R cos 0 R g> R - sin $ R 0 r 


[GB] ^ 
r I 


L^ tan ♦cj 
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Integration Variables 

The number of integrals computed during any particular phase 
is determined from the options requested by the user. As a mln- 
imum, the translational equations of motion are Integrated to 
give the position and velocity of the center of mass of the ve-= 
hide. The user may also select additional vat tables to be in- 
tegrated. The only restriction is that no more than AO integrals 
can be computed per phase. 



! 
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VII. AUXILIARY CALCULATIONS 


In addition to computing the basic variables, POST also com- 
putes numerous auxiliary variables that are related to: (1) 

conic parameters, (2) range calculations, (3) tracking data, (4) 
analytic impact calculations, (5) velocity losses, and (6) veloc- 
ity margins. The equations used to calculate these variables are 
presented below. 


Conic Calculations 

The following Keplerian conic variables are computed. 

♦ 

V 2 

I ji_ 

energy, ^ 


semimajor axis, -u/26 
angular momentum, |xj x J2jl 

semilatus rectum, h 2 /v 
eccentricity, /|1 - p/a| 

velocity required to circularize orbit, /AV • AV, where 
u^ ■ h/h 

%i ■ £ i /r i 


Hv -Sh '.Stl'ISh * Sri 1 

V. - (M/r,) 15 !!, 


4V ■ V - V, 


inclination, cos" 1 ( h 8 /h) 

longitude of ascending node, cos* 1 (Xj • u^), where 


VII-1 


■ z. j * h / IjSj * h 


argument of vehicle, p » cos" 1 ^ u^, • u^ 


time since perigee, ^ M 


time to perigee, P - T 


SP 


latitude of perigee, tan -1 (03/ /ij + uf) » where 
u * eosfw)^ + ain(o))|u^ * u ( ^ 

longitude of perigee, tan" 1 (U2/U1) 


altitude of perigee, r - R ($ ) 

P 8 P 

altitude of apogee, r - R ($ ) 

asp 


velocity at perigee, 


>. J \ (Hi) 

velocity at apogee, J ^ 
hyperbolic excess velocity, / 26 

maximum true anomaly for hyperbolic orbit, cos" 1 (-1/e) 
declination of outgoing asymptote, sin" 1 lu y ( 3 ) 1 , where 

u) 

St ■ Sh “ Sri 

»r ' COB(9 m« - e) Sri + » tn < 8 MJI - 9 > Sr 


right ascension of outgoing asymptote 




I 


f 


i 


i 


M 

0 ) 

r 


true anomaly, cos" 1 - ljj 

eccentric anomaly, 2 


tan 


!) 


mean anomaly, & - e sin E 
argument of perigee, p - 0 
perigee radius, a(l - e) 
apogee radius, a(l + e) 


period 




% 
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Range Calculations 

The progam provides for various types of range calculations. 
The equations for these calculations are given below. 

Dot product downrange .- The relative range angle, measured 
from the vehicle's Initial position to its current position, is 
given by 


* 


R 



(VII-l) 


where u^ is a unit vector along the initial position vector 
so 

in Earth-centered rotating coordinates and u is a unit vector 

along the current position vector in Earth-centered rotating co- 
ordinates. The range over an oblate spheroid is calculated from 
the average radius to the surface, and is given by 


*D 



(VII-2) 


Crossrange and downrange via orbital plane reference .- Re- 
f erring to figure VII-1, identify the vehicle's position at time 

^ A 

t by 0, and at a later time t by P. At time t , the ve- 
hicle has a latitude of $ , a longitude of 0 , and a velocity 
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F~ r ~ - r 



o 


<0 


0 
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heading of X . At time t the vehicle Is at latitude $ and 
longitude 6. The downrange angle (y) and the crosstange angle 
(v) shown in the illustration are measured along, and normal to, 
the great circle through 0, and are inclined to the meridian by 
X*. From analytical geometry, v and u can be expressed as 

^ jUl jUi 

sin v ■ -sin X sin $ cos $ cos 0' - cos X cos 0 sin 0 

c c 

* * 

+ sin X cos $ sin 

w 

sin y ■ (-cos X sin 4> cos $ cos 0' + sin X cos <J> c sin 0' 

+ cos X cos $ sin $ c j/cos v 
cos y “ (cos 6 cos 4 cos 0" + sin $ sin 4> c j/cos v, . 

ft ft 

where 0 and X can be defined in either of two ways: 

1) The great circle to which v and y are referenced Is 
fixed and rotating with the Earth. Then 
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X* « Earth's relative heading » sin” 1 


r - 8 - 0 ; 


J4 


+ ti 


? (V1I-4) 


2) The great circle to which y and v are referenced is 
inertlally fixed, having the Earth rotating below If. 
Then 


inertial heading - sin -1 


0' - 0 - 0* + fl (t - t*). 
P 


+ v 2 


(VII-5) 
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Auxiliary Position and Velocity Calculations 


The solution from the translational equations is then used 
to calculate numerous output variables. The key variables directly 
computed from (Xj* y^» Zj) and (V XJ ., V^j, V ZI ) summarized below. 


h 


h 


'R 


S*I 


% I 


geocentric radius 
*3 


(* • *) 
magnitude of the inertial velocity 

(2l • li) H 

relative velocity 

V «• 0 x r 
-I Zf * -I 

atmospheric relative velocity 


V + V 


magnitude of the relative velocity 

(s* • h? 


\ 


magnitude of the atmospheric relative velocity 

(4 • s*) % 

unit vector along radius vector 

hl'i 

unit vector along inertial velocity vector 
i 

inertial flight path angle 
sin“ l 


[s*i - *vt] 


> (VII- 


■ relative-flight path angle 

' “°" 1 [a«i • V] 

Y a - atmospheric relative flight path anjjle 

- •*»-* [%, • Uv A ] 

2io * Inertial velocity in the G-frarae 

- [IGJ Vj. 

^RG relative velocity in the G-fram6 

- [IG] V* 

-AG atmospheric relative velocity in the G-frame 

* tIG] V A 

A ZI Inertial azimuth 

- «“•' [ v yg/ v xc] 

A^ » relative azimuth 

■ [ v ryg/ v rxc] 

* atmospheric relative azimuth 

■ t “' 1 [ v ayg/ v axg] 

♦ - geocentric latitude 

* * ln ' 1 [*l/ r t] 

0j ■ inertial longlt 'e 

■ t,n ‘‘ [Vi] 
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e R ■ relative longitude 

■ "i - S (' - e o) 

Agg * sensed acceleration In the B-frame . 

^TB + — AB 

Ag ■ magnitude of the sensed acceleration 

"(is • is)* 

Agj ■ sensed acceleration in the ECl-frame 

' '“1“ [4r> + 4 aj] 


> (VIX-7) 
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Auxilary Attitude Calculations 


The attitude angles that are not used to generate the steer- 
ing commands are computed for output In the auxiliary calculation 
subroutine. These equations are summarized below: 

1) Aerodynamic angles: 


a 

6 

a 


tan' 


-1 


tan' 


-1 


tan 


-l 


(* B / U B ) 

( V B^/ U B ♦ W B 2 )* 

( GB 2 3 + sin $ sin^ 

GB 22 co* - GB 2 i sin A 2A cos 


> (VII-8) 


2) Inertial Euler angles: 


♦i “ 

tan -1 

(lb 23 /lb 22 ) , 

0 

M 

-sin -1 

(lb 2 i). 

9 i a 

tan” 1 

(lb 31 /lb h ); 

3) Relative 

Euler angles: 

'“r - 

tan -1 

(GB 12 /fcB u ), 

6 R- 

-sin" 

1 (gb 13 ). 

♦r" 

tan" 1 

( GB 23/ GB 33) • 


1 OTI-9) 
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Tracking Data 


?08T computes tracking Information for as many as tan 
tracking stations per phase. The tracking stations are located 
on a reference ellipsoid and are specified In terms of their lati- 
tude* longitude* and altitude above the ellipsoid. These variables 
are Illustrated In figure VX1-2. 



Figure VI1-2.- Radar Tracking Schematic 
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The position components of the tracker in the Earth-relative 
frame are given by 


%R 



*cos <|> T cos t» T * 
cos <> T sin » T , 


(VII-11) 


where ttj. Is the altitude of the tracker, ^ the latitude of 
the tracker, and the longitude of the tracker. 


The slant range vector in the ECI frame is given by 

£ sr - Ij - UP)- 1 It R . (VI1-12) 

and the slant range is then computed as 


' S R • -J '- SR • £SR- 

The elevation angle can then be computed as 

’t • (% • Jsr/^r)' (VII-W 

where 

iijR - [IP]- 1 L m / I [IP]” 1 £ TR | . (Vll-15) 

The slant range vector, transformed to the geographic frame, 
is 

XSRC * ‘“1 isR’ < VII -“> 

and thus the tracker's azimuth is given by 

*2T ■ ( y SR0/ X SRc)' <VI1 - 17 > 

The look angles are calculated from the slant range vector 
transformed to the body frame; l.e., 

Isrb - [IB] rg R . (VIl-18, 
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~T 


0 
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Using the components of £g RB * the eone angle ia than given by 

- cos - 1 (* SRB / r SR )» 

and the clock angle Is given by 

“c • ton ‘‘ ( ^SRB /*SRB ) ‘ (VIl-20) 

Space losses are calculated for the tracking stations as 
follows: 

SLj - 36.56 + 20 Logic (R gLM • F R i) 

SLj * 36.56 + 20 Logjo ^ R g y>j * FRz) 

SL 3 - 36.56 + 20 Log io ( R SLM * FR3) . 

where 

FRi ■ 420.0 (command frequency) 

fr 2 ■ 2287.5 (telemetry frequency) 

FR 3 * 5765.0 (tracking frequency) 

R_ „ ■ slant range distance In statute miles. 

SLM 

Analytic Impact Calculations 

The analytic impact calculations predict the geodetic lati- 
tude, longitude, and time of flight at impact for a vehicle with 
a given position and velocity to its intersection with the sur- 
face of the oblate planet. These calculations assume Keplerlan 
motion and are not corrected for drag effects. 

The basic problem in determining an impact point from a 
specified position and velocity |r IQ , V. I0 ) is in calculating 

the Impact eccentric anomaly. This angle is determined by 
iteratively solving the equation 


VII-13 



(VII-19) 




+ ll 


ip 


where h 

ip 

planet and 


In the denlred impact altitude above the oblate 
tin* partition vector l» given by 


(VII-22) 


,<K) « 

•'<« - 1 y l0 


1(F) " 

(eon (K - I'.q) -o cos Kq) j (l - e coh Ky) 

(VII-23) 

R(K) 

/""p- (a In E ~ Eyj -e sin 1. + e sin Ky . 



the 


Once the impact eccentric anomaly, E^, * s determined, then 
time, latitude, and longitude of impact arc calculated as 



(VI 1-2 A) 
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VIII. TARGETING- AND OPTIMIZATION 


POST uses an accelerated projected gradient algorithm (PGA) 
ae the baalc targetlng/optlmlzatlon technique. PGA le a combina- 
tion of Rosen’ a projection method for nonlinear programming (ref a. 
3, 4, and 5) and David on's variable metric method for unconstrained 
optimization (ref . 6) . The program also contains backup slngle- 
ponalty function methods that use steepest descent , conjugate 
gradients, and/or the Davldon method. These standard gradient 
methods are well documented In references 6 and 7 and are only 
briefly described in the following discussion. 

The projected gradient algorithm is an Iterative technique 
designed to solve a general Class of nonlinear programming prob- 
lems. PGA employs cost-function end constraint gradient Informa- 
tion to replace the multidimensional optimization problem by an 
equivalent sequence of one-dimensional searches. In this manner. 

It solves a difficult multidimensional problem by solving a se- 
quence of simpler problems. In general, at the Initiation of the 
Iteration sequence, PGA Is primarily a constralnt-satisflcation 
algorithm. As the iteration process proceeds, the emphasis 
changes from constraint satisfaction to cost-function reduction. 

The logic used to effect this changeover process will be dis- 
cussed below. 

Since numerous analytical developments of this technique are 
available (see refs. 3, 4, and 3), this presentation will pri- 
marily emphasize the geometrical aspects of the algorithm. This 
geometric Interpretation clearly motivates the equations and 
logic contained in PGA, and a basic understanding of these con- 
cepts is usually sufficient to enable the user to efficiently 
use the algorithm. 


Problem Formulation 

The projected gradient method solves the following nonlinear 
programming problem: 

Determine the values of the Independent variables, u, that mini- 
mize the cost function (optimization variable) 

F(u) , (VIII-1) 


! 
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* r 




o 






(VIII-2) 


subject to the constraints (dependent variables) 

c(u) 0, 

where u c j c is a vector - valued function* i.e., _c»R R * 
and F is a scalar-valued function, i.e., F:R * R . 

The algorithm is actually more versatile than this simple 
formulation might indicate. In order to maximize any particular 
function, say W(u), all that is required is to define 
F(u) * -W(u) and determine the minimum of F(u). The equality 
constraint - case is also contained within the above formulation 
since constraint equations of the form 

c (u) - 0 (VIII-3) 

J 

are special cases of Eq (VI11-2). 

In the trajectory optimization, the cost function end the 
constraints are not explicitly a function of the independent 
variables, but rather depend explicitly on the stae variables 
r Vj, m, and Q. The explicity equations relating the state 

(dependent) variables to the independent variables are the in- 
tegrals 


- r -I 




dt 


l z “ 1 1q + + 4b] + £l] dt 



(VIII-4) 


If x denotes the above etate varlablea-of the eyetem being 
alaulated at the n event, and x and x denote the value 

"ti — n 

of on the plus and minus sides of that event, then 

Vi ‘ sj . mn-n 

where are the Independent variables lr. phase n, and 

represent the solution of the state differential equations over 
phase n. The values of the state variables on the positive side 
of event n are them 


+ 

*n 


+ Ax 
~n 


(vin-e) 


where Ax represents the discontinuity in state (e.g. , velocity 
“ th 

Impulse at the n event) . 

The cost function and the trajectory constraints are computed 
at the positive side of the specified events, and are therefore 
given by 

(VIII-7) 


(VIII-8) 


la specified and v^ denotes the events at which the dependent 

variables are specified. This generality enables the program to 
solve problems in which intermediate constraints are defined, as 
well as problems where the cost function is not specified at the 
final event. 


F(u) 


u> - <(*;,) 


and 


I “V. 


i K) 


C (u) - 


l-V 


3 Kj) 


where v. denotes the event at which tha optimisation variable 
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The trajectory propagator, T n » can represent either numer- 
ical Integration or analytical Keplerian equations. 


Fundamental Concepts and Homendature 

To facilitate the discussion of the projected gradient algo- 
rithm, the following nomenclature and basic concepts will be In- 
troduced . 

A real k-dimensional Euclidean vector space Is denoted by 

R k , and x denotes a column matrix whose elements ate 
x , where 1 » 1, 2, k. The vector Inequality x >_ 0 Im- 
plies x^ 0 for each i, and A* denotes the transpose of the 

real matrix A. 

The cost gradient Is an m-vector of partial derivatives de- 
noted as VF or 3F/9u, and Is defined as 


^ ■ FT 


(VIII-9) 


The gradient to the 1 th constraint Is similarly represented. 


The Jacobian matrix of the constraint vector function with 

respect to the independent variable Is a matrix whose 1 th row Is 
the gradient vector Vc^. This matrix is denoted as 

9c 

J(u) - r- (VIII-10) 


and contains n rows and m columns. Clearly, 


ij 8u. 


(VIII-11) 
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The J^ h constraint is said- to. be active at u if and only if 

a) c^n) < 0, (VI 11-12) 

An active constraint is said to be unconstraining if and only if 

b) Cj (u) - 0 and » [(SS') _1 1 0. (VIII-13) 

Condition a) implies that the j ” constraint is either vio- 
lated at u, while b) indicates that the negative of the cost 
function gradient "points" outside the feasible region. 

The sensitivity matrix is that matrix whose rows are the 
gradients to the active constraints, and is denoted by 

3e 

S(u) - (VIII-U) 

where e is the n ft -vector of active constraints. Equality con- 
straints are always active and thus are loaded into the upper 
elements of the e. Thus, e is essentially the error vector 
for the active constraints. The error function is defined to be 

E(u) - e'e. (VIII-15) 

The sensitivity matrix, S, is obtained from the Jacobian 
matrix, J, simply by deleting those rows that correspond to 
Inactive constraints. 

Corresponding to each constraint function c^u) is a 
boundary hypersurfaoe , B^» defined by 

B i " jH sc i<H.) • oj. (VIII-16) 

Clearly, B 1 Is an m-1 dimensional nonlinear manifold. It can, 

however, be approximated at each point 4 in R 11 by an m-1 
dimensional linear manifold 

C^(d)* (u - d) + c^(u) * o|. (VXII-17) 
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The feasible region for the i*** Inequality constraint la the 
half-space In the Independent-variable space defined by the set 

\ 1 0}. (VIII-18) 

while the complete feasible region for all of the constraints Is 

n 

R- O R a . (VIII-19) 

i“l . 

The boundary of the complete feasible region must be 

n 

MR) * U (b t n R) (VIII-20) 

1-1 

The Intersection In the preceding equation Is required .to select 
from the unbounded boundary, B^ of the feasible region of the 

i** 1 constraint that portion which is adjacent to the feasible 
region, R, for all of the constraints. 

At any particular CL k R m it is useful to define the local 
boundary hypersurface, B(A), .to the complete feasible region as 
the intersection of the active constraints at 4. Let N(tt) 
denote the set of indices of the tight constraints at A. 

Then, symbolically, 

B(A) “ O (VIII-21) 

ieK(u) 

Clearly B(A) is an m-R dimensional nonlinear manifold in the 

a 

m-dimenslonal independent variable space* 

An m-k dimensional linear manifold C(A) approximating 
a 

B(fl) is the intersection of the active linearised constraints at 
A; that is. 
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o 


(VXII-22) 


c(0) - n ^(u) 

lcK(O) 


- |u:S(0) (u - d) + e(d) - o| 


(VIII-23) 


Mow let ?Ka) denote the llneer spece spanned by the gradients 
to the active constraints; that Is, 


k a l 

$<u) u:3 a x a n for which 

l a 1-1 

and let Q(d) denote the orthogonal complement to §(<*)» that 

Is, 


(VIII-24) 


R m - Q(d) © $(d). (VIII-25) 

It can be shown that Q(d) Is the unique linear space that can 
be translated to obtain the linear manifold C(d) . 

Furthermore there exist unique orthogonal projection oper- 
ators P(0) and ?(u) that resolve any vector in the independent- 
variable space Into Its corresponding components in Q(d) and 


$<i). 

respectively; tha.t Is, for any u c R m 



u - p (u)u + $(d)u» 

(VI 11-26) 

where 

P(d)u e Q(d) and l(0)u e $(u) . 

(VIII-27) 


In particular, 



- S^SS-)" 1 3 

(VIII-28) 

and 

P - I - 

(VIII-29) 


An additional concept Is the Idea of problem scaling. The 
purpose of problem scaling la to Increase the efficiency of the 
targeting/optimisation algorithms by transforming the original 
problem into an equivalent problem that is numerically easier to 

solve. 

To numerically scale a problem, two general types of scaling 
are required: (1) independent-variable scaling, and (2) dependent- 

varlablu scaling. Independent-variable scaling is accomplished 
by defining a po- tive diagonal scaling matrix, W^, such that, 

the weighted inJepend nt vcxi'vab lee are given by 



(VIIT-30) 


Simularly, dependent-variable weighting is accomplished by 
defining an optimization-variable scale factor, W f » an 

positive, diagonal, dependent-variable scaling matrix, such 

that the weighted optimization variable is 

p . W _ F (u) (VIII-31) 

X F "" 


the weighted dependent variables are given by 


£<£> ■ [«.j£ (V‘ a )’ 

yielding a weighted error function 

p 2 - • • (y) ♦ 


(VIII-32) 


(VIII-33) 


The program contains several options for computing the in- 
dependent-variable weighting matrix. However, the option most 
often used is the percentage scaling matrix 



(VIII-34) 
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The dependent -variable weighting matrix la always computed 
as the reciprocal of the constraint tolerances, and Is given, by 



(VIII-35) 


where t ± Is the tolerance for the 1 th constraint* The optimisa- 
tion scale factor Is merely Input so that Pj Is approximately 
equal to one. 


For simplicity, the following discussion of the algorithm 
assumes an appropriately scaled problem* However, the scaled 
equations can be obtained by making the following simple sub- 
stitutions: 


u replaced by ti 
F replaced by P^ 

x 

c replaced by c 


h £ replaced by Pj 
S replaced by JwJ [S] |w^ * 
VF replaced by W (| 1 £F* 


The final key concept employed by PGA is the Idea of a direc- 
tion of search. Heuristlcally, the direction of search is nothing 
more than a particular line In the independent-variable space 
along which the constraint error Is reduced, or along which the 
cost-function Is decreased. In a more precise sense, the direc- 
tion- of search at fl Is a half -ray emanating from fl. Thus, for 
any positive scalar, y, the equation 

u - Q + y! (VIII-36) 

sets the limits of this half-ray and represents "movement" in the 
direction i from 0. This Is Illustrated In figure VIII -1. 
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Figute VIII-1.- Direction of Search in the 
Independent-Variable Space 

If £ is a unit vector, then •» represents the actual dis- 
tance "moved " in the direction This concept of direction-of- 

search is particular 1.. important since it enables the m-dimen- 
sional nonlinear programming problem to be replaced by a sequence 
(hopefully finite) of one-dimensional minimizations. What remains 
to be explained then is: (1) how to select the direction-of- 

search; and (2) how to determine the step size in that direction. 
All "direct" optimization methods employ this concept and, hence, 
differ only in their answers to the two preceding questions. The 
technique by which and > n are selected by PGA will be de- 

scribed in subsequent sections. 


Direction of Search 

The projected gradient method uses two basic search direc- 
tions. For this discussion, they will be termed the constraint 
and optimization directions, respectively. PGA proceeds. by tak- 
ing successive steps in one or the other of these two directions. 
The computation of each of these search directions is described 
below at a particular point d in the independent-variable space 
where A of the constraints are active. 


V1II-10 
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Constraint d tract ion .- Tha constraint direction depends 
critically on the number of active constraints* Three cases ere 
distinguished bslowj 

1) Case 1.- If m, then that unique c* .trol correc- 

tion AO is sought, which solves the linearised con- 
straint equation 

S(ft) Au + e(d) - 0 (VIII-37) 

and minlalzes the length of Au. The solutions to the 
preceding vector equations defTne the m-ft dimensional 

A 

linear manifold C(fl) , which approximates the local 
boundary at £ as described in detail in the preceding 
section. The desired minimum norm correction, Ad, is 
then the vector of minimum length in the indepdendent- 
vavlable space from d to the linear manifold C(d). 
Analytically, it is given as 


AO - -S '[ SS ']"' 1 e(d) . (VIII-38) 

This correction Is Illustrated in figure VIII-2. 

The direction of search then la simply taken to be this 
minimum-norm correction to the locally active linearised 
constraints; that is, 

s c (d) - A$. (VXII-39) 



Figure VIII-2*- Illustration of Minimum-Norm Constraint, 
Direction for ■ 2 < m ■ 3 
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If d m m , than the linearised local boundary C(a) 

a 

reduces Co a single point. Thus, chare Is a unique solu- 
tion to the linearized constraint equations without the 
additional requirement chat the length of the Independent- 
variable correction be minimized. The minimum-norm cor- 
rection formula then reduces to the familar Newton- 
Raphdon formula for solving m equations In m unknowns; 
nameJ.y 

Ah - - S" 1 e(Q) . (VII] 

The Newton-Raphson correction Is Illustrated geometrically 
in figure V11I-3. 


Second linearized 
constraint 


Third linearized 
constraint 


Ad, Newton-Raphson 
correction 



First linearized 
constraint 



c(d), intersection of 
linearized constraints 


Figure VIII-3.- Illustration of Newton-Raphson Constraint, Direction 


for n. - m - 3 


The direction of search is taken to be this unique cor- 
rection vector satisfying the linearized constraints; 
that is, 


s(d) - Ad. 


(Vlil-41' 
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(V1II-43) 
>n pic 









Optimization direction .- When the number of active conetralnte 
ie leaf* than the number of Independent variables, it is then pos- 
sible to reduce the nonminimal cost-function. Obviously the 
steepest descent direction, - VF (il ) , would be the best local 
search direction for reducing the cost function. Such a direc- 
tion, however, would generally produce unacceptable constraint 
violations. To avoid this difficulty PGA orthogonally projects 
the unconstrained negative gradient, -V F (0 ) , into a direction 
parallel to the local linearized constraint boundary C(d>. By 
searching in the direction of this negative-projected gradient 
the algorithm can guarantee that there is no further constraint 
violation than that of d for the case of linear constraints. 

To calculate this direction, it is only necessary to apply to 
the unconstrained negative gradient the projection operator P(u)» 
which maps any vector in the independent -variable space into its 
component in Q(d) , the unique linear space that can be trans- 
lated into coincidence with the linear manifold C<d) . Thus, 

s°(d) = -P(u) VF(d) \ 

« -[I - P(d) ] VF(d) HVIII-45) 

- -[I -S' (SSO“ l S(d)l VF(d) J 

The direction of search for the accelerated projected gradient 
method is 

s°(d) * -H P VF(fl) (VIII-46) 

~n n — 


where 

H o * 1 


H »- 

Vi + \ ♦ V 

where n ■ 2 

A n‘ 

b* **]/** 

v 

®n " 



A *n“ 

Wi’ 


*n“ 


~l)‘ 


(VIII-47) 


(VIII-48) 
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figure VIII-5 illustrates the direction of the negative-projected 
gradient for the case-o£-a_elngle active constraint. 


first ioMttvi 
constraint 


Lines r Usd approximation 
to sols set i vs constraint 

H * 

Unconstrained negative 
gradient to cost 
function st £ 


Msgs ti vs projected 
gradient st w 



Sole setive , 
constraint 
(paraboloid 
of revolution) 


Second inactive 
constraint 


Figure VIII-5.*- Direction of Negative-Projected Gradient for n ft - 1 

a nd m « 3 (Feasible region is that region inside 
paraboloid, above lower plane, and below upper plane; 
coet-functlon is vertical height) 

if there are no equality conatraints, and if all the inequality 
constraints are Inactive, then S is the zero matrix and the 
direction of search becomes the standard deflected gradient 
direction 


s° (fi) - -H n 7F<a). (VIII-49) 

Similarly, if the single-pefcalty-f unction methods are used, 
then the directions of search that minimize 

P 2 - F + W e 'e (VIII-50) 

are; 


1) Steepest-descent method 
s°(fl) - -VP 2 (fl); 
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2) Conjugate gradient method (steepfiAfcrd.e.8pent starter) 

s , , where n > 2 ; 
— n-1 ' 


s° - -VP, u + 
— n — 2 — n 


VP 

2 1 

r” 

)*2 



■x) 


Vi 


3) Davidon's method (steepest-descent starter) 

s° - -H VP,(u \ , where n >. 2 
— n n — 2\-n/ 


and 


H ■ H . + A + B , 
n n-1 n n 


where A and B have the same definitions as in the 
n n 

accelerated projected gradient mode. 


Step-Size Calculation 

At any particular point £ in the independent-variable 
space, the PGA algorithm proceeds by reducing the multidimen- 
8ional problem to a one-dimensional search along the constraint 
direction to minimize the sum of the- squares of the constraint 
violations, or along the optimization direction to minimize the 
estimated net cost-function. In either case, once the Initial 
point d and the direction of search £ are specified, the prob- 
lem reduces to the numerical minimization of a function of a single 
variable— namely, the step size. PGA performs this numerical 
minimization via polynominal interpolation, based on function 
values along the search ray and the function’s value and slope at 
the starting point. Consider then, in detail, the calculation o 
this latter pair of quantities for the respective functions asso- 
ciated with the constraint and optimization directions. 


Constraint direction . - The function to be minimized along 

the constraint direction, £ c , is the sum of the squares of the 
constraint violations; namely 

h c (Y> - li(a+vi c )l 2 - (viii-sn 
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Clearly 

h (0) - |e(0)| 2 . (V1II-52) 

c 

Differentiation via the chain rule yields 

h c '<0) - 2e'(a)S(a)i C . (VIII-53) 

Recall that the search direction s c was obtained as an in- 
dependent-variable correction either satisfying all the linearized 
constraint equations If n^ £ n, or minimizing their violation if 

m < n . Thus, If the constraints are reasonably linear, a good 
a 

initial estimate for the y minimizing h, is one. 


Optimization direction .- The function to be minimized along 

the optimization direction, s°. Is the estimated net cost- 
functiort which Is defined ad 


h o (y) ■ f(o + ye°) - *(<*) 



change in cost- 
function produced . 
by step of length 

y along 8° 



linearized approximation to 
change in cost-function re- 
quired to perform minimum- 
norm correction back to the 
feasible region 


(VIII-54) 


Clearly 

h o (0) - -r?(<9S'(SS'r l <£)«(£)• (VIII-55) 

By expanding h Q in a Taylor aeries In y about y - 0, 

and by making use of the fact that $a° - 0 since i° lies In 
Q(fl) , it can be shown that 

h'<0) - s°. (VIII-56) 

These properties of h Q are Illustrated In figure 25. 
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Figure VIII-6.- Properties of Estimated Net Cost Function 

Both the constraint and optimization directions are based on 
a sensitivity matrix that depends critically on which constraints 
are active. Hence, for searches in either direction, it is im- 
portant to limit the step size so that the set of active constraints 
does not grow. Such a limit can be obtained based on linear ap- 
proximation and suffices to deal wit' inactive constraints becom- 
ing active. 

The reverse situation — of active constraints becoming in- 
active— poses no difficulty. To see this, note that because of 
our treatment of the active constraints as linear manifolds, a 
first-order approximation of the distance to a particular active 
constraint boundary would not change along the optimization direc- 
tion. Furthermore, along the constraint direction any change in 
the status of an active constraint will be appropriately treated 
by minimizing h with respect to the step length. 
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Lit K(fl) dmota thi ait of activa constraint Indices at 
4, and let 

r k - a'(4)Vc k <4) , (VIII-57) 


where £(4) Is the search direction at vector 4. Then assign 
to each” k In K the number 


A(k) 


-c k (4> /c k if t k < o 

R If r k > 0 


I 


(VIII-58) 


where R Is a very large real number* Then A(k) Is a linear 
approximation to the distance along the search ray from u to 

the boundary, B k , of the k ttl cotostralnt. Hence a resonable 

upper bound for the step length Is 

X - min [A(k)J * (V1II-59) 

kcK 


One-Dimensional Minimisation 

Monovariant minimisation In PGA is performed exclusively by 
polynomial Interpolation. First the actual function, f, to be 
i» fitted with one or more quadratic or cubic poly- 
nomials until a sufficiently accurate curve fit, p, is ob- 
tained; that is, 

n 

p( Y ) a 1 Y i «^f(Y) for all y of Interest. (V1II-60) 

1-0 

Then the Independent variable value, y B , that minimises f Is 

approximated by the value, y B » which minimises p. Clearly, 

y® can be- determined analytically If n <. 3. 

P 
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The minimization routine makes ingenious use of all the in- 
formation it accumulates about f to obtain a good curve fit. 
First, f is fitted with a quadratic polynominal, Pi» based 
on: 

1 ) *< 0 ) 

2) f'(0) 

3) where y* > 0 is an initial estimate of the 
Y value that minimizes f. 


The coefficients of this quadratic polynominal are then calcu- 
lated from the formulas: 


The value of 
nominal is 


ao • f(0) 
aj * f '(0) 

12 ‘ [ r ( , o) • *«]/ T ' 2 + *'/ '“• 

the independent variable that minimizes 


| (VIII-U.) 
this poly- 


Y ® - -a x /2a 2 . (VIII-62) 

If y? and Yo do not differ significantly, y”* is taken 
to be y? and the minimization procedure is considered complete. 
Similarly, if Pi (>?) * 8 not significantly different from 

f (yT)» then is taken to be equal to y? and the process 
is terminated. Otherwise f is fitted with a cubic polynominal, 

P 2 » based on 

1) f(0) 

2) f'(0) 

3) f (yo) and Vo > 0 

4) f (y?). 
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If f la fitted using P 2 » then coefficients ere calculated 
from the following formulas: 

*\ 

a 0 — f(0) 

Si - f'(0) 

X - max (y?,vT) 
a - min (yo»Y?)^ 

*2 ■ [Xaj a + ao (1 + ot) + (a 2 f(X) “ f(«X))/(l - ot)]/(X 2 a 2 ) 

*3 * Kf(aX) - <* 3 f(X»/(l - «)-Xa(l t a)ai-(l + o - a 2 )a 0 ]/(X 2 « 2 )^ 

The value of the Independent variable, X™,- that minimises this 
cubic polynomial is 

Y® * (“®2 ♦ \^2 “ 3 * 3 * 1 ) / 3 *3 • 

If Y 2 and Y? do not differ significantly, y® is taken 
to be Y2 40(1 the minimisation Is stopped. Similarly, if P2 (y® ) 
is not significantly different from f (y®) » then y® 1® taken 
to be equal to y® and the procedure is terminated. 

If none of these stopping conditions Is met, a third quad- 
ratic curve-fit Is attempted. The accumulated set of sample 

nointa on “ f , namely [0,f(0>], £y®, f (y®) [Y?,f (y?)J . and 

Ty 2 » f (y?) 1 , Is arranged In the order of their ascendiiig abscissa 
values . Then the first point whose ordinate value is less than 
that of the following point Is selected. 

To simplify the notation in the following pages., relable this 
point as Iy 2» f(Y2^1» the preceding point as [yi» f(Yi)l» 
the following point as [ys» f(Ys)]* 


(VIII-63) 
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Another quadratic polynomial, 


la then fitted to 


P3 » 


1) 

f<Yl) 

2) 

f(Y 2 ) 

3) 

f(Y3>‘ 


The formulas for these quadratic coefficients are as follows: 



" Vj 






c 13 

* Y i + Yj 







“ Y i - h 







b 2 3 

f(u) 

bi3 

f(Y 2 ) 

b 12 

f(Y3> 

ao 

di 2 di3 

d 2 id 2 3 

d 3 1 d 3 2 


c 23 

f(Yi) 

C 1 3 

f (y 2 ) 

c 12 

f(Y3> 

•l 

di 2 dl3 

d 2 id 2 3 

d3id 32 


1 

f.(y i ) 

+ 1 

f(Y 2 ) 

+ 1 

i(Y3>‘ 

a 2 

dj 2 di3 

d 2 id 2 3 

d31 d 32 
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■ The value of the independent variable that minimises this quad' 
ratic la 


Y ? - - ai /2a 2 . (VII 1-66) 

If y” and y“ do not differ significantly, Y m is taken 
to be y? and the search ip discontinued. On the other hand, if 
p 3 (y “ ) is not significantly different from f \Y 3 ) » then y 
Is taken to be (y®) and the process is terminated. 

If neither of these stopping conditions is met, then a cubic 
polynomial is fitted to 
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D f(Yl). Yl 

2 ) f(Y2>. Y 2 

3 ) f(Y3>. Y 3 

4 ) f(Yi*)» Y«4 " 


The formulas for these coefficients are as follows: 
Di - (Y 2 - YlHrs - Y 1 ) (Y 4 - Yl) 

D 2 ■ (Yl - Y 2 ) (Y 3 " Y2HY4 - Y 2 ) 

D 3 • (Yl - Y3XY2 “ Y3HY4 " Y3) 


D4 • (Yl - Y4) (Y 3 ~ >4)(Y 3 “ Y4) 


Y2Y3Y4 


Y1V3Y4 


Y 1 Y 2 Y 4 


Y 1 > 2 Y 3 


*0 “ pf - f(Y ^ + ' D 2 ~ f < Y 2> + D 3 ” f(Y3> + D4 f(V4) 

Y 2 Y 3 + Y 2 Y4 + Y3Y4 (YlY 3 + Yl *4 + Y4Y3) 


Y2Y3 + Y 2 Y 4 + Y3Y4 
*1 x £(yi) ♦ 


£(y 2 ) 


(YiY 2 + Y1Y4 + Y 2 Y 4 ^ (Y1Y2 + Y lY 3 + Y2Y3) 

+ 5 f(Y 3 > + Dj f(Y4) 

(Y 2 + Y 3 ♦ Y43 (Yl + Y3 + Y4) 

•2 - i r 1 f(Yl) + 5:; f <Y 2 ) 


(Yl + Y2 + Y 4 > 


f(Y 3 ) + 


(ri + Y2 + Y4) 


f(Y 4 > 


*3 * " dT f<Yl> " f<V2) " »| £<Y3) “ f< Y4 ) * J 

The value of tho indp undent variable minimizing this fourth cubic 
polynomial la 

Y4 • (-*2 t / a 2 - 3 a 3 ai)/ 3 a 3 . 
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If y? and y! do not differ significantly, y n la taken 
to be y® and the minimization ls-stopped. Similarly, If 
Pi^Y®) 1* not significantly different from f (y™) • then y n 
Is taken to be equal to y? end the procedure Is terminated, 

If none of these stopplug conditions Is met, the accumulated 
set of sample points is searched for the point with the minimum 
ordinate value. The abscissa value of this point Is taken to be 

Y°, and the minimization Is considered complete. 


Algorithm Macrologic 

After being Initialized the projected gradient algorithm 
proceeds as a sequence of Iterations, each consisting of an op- 
timization step followed by a constraint-correction step (sea 
fig. VIII-7) . The very tirst step from the user's initial independ- 
ent-variable estimate is however, one of constraint correction. 
Furthermore, the optimization step Is also omitted on any itera- 
tion for which the constraint-violation function, h , was not 

c 

reduced by the constraint correction step of the preceding Itera- 
tion. 


The optimization search direction that emanates for u^ is 
based on the sensitivity matrix, S fu \ ; that Is, 


4 * a °(^n) ' - p * F (M • 

as discussed previously. Hence, s£ lies in the 

The value of the Independent -variable vector, 
the optimization le 


(VIII-69) 
subspace Q^u^j . 
u° , after 


o o 

v , • JL + 
f\ *11 o — n 

where y 0 is the optimum step site. 
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sr - 
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Flgur« VIII-7.- 




Macro logic of Projected 
Gradient Algorithm 
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I 

I 

1 


The direction of the constraint-correction search emanates 
from u°; however, since generating a new senaltlvlty matrix 

—tl 

is such an expensive calculation, the old Jacobian matrix, J, 
of the constraints with respect to the controls evaluated at 

u Is used In conjunction with the error at u°. Thus, 


^ - -S'(SS-)-*(u n )s(u°) 


It can be shown by direct computation that 


*(u\s c - s< 
y-n)-n -r 


where 


space 


Is based -on S(u Thus, s c lies in the sub- 
\—n/ ’ — n 

in the independent -variable space. 


Since and Q^u^j are orthogonal complements. It 

follows that the optimization and constraint directions for any 
iteration are exactly orthogonal; that Is, 


to' 


s c » 0. 
— n 


The result of the constraint correction step Is then the Inde- 
pendent-variable vector for the next Iteration. Thus 

*n + i - 4 + Y c V 


Figure V1II-8 


PGA Iteration. 
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Unconstrained gradient 
of coot function, £F(u) 


Plane determined by 
gradient to the coat 
function and the gradient 
to the active constraint 

/— Minimum-norm constraint 


Gradient to c^, , 
Vcj (G) 


Linearization of 
sole active 
constraint 


2 $ 


Sola active nonlinear 
constraint, C 


Projected gradient { 
optimization step, y s 


Figure VIII-8.- Complete PGA Iteration, Consisting of Optimization 
* Step Followed by Constraint Step for n fl » 1 and 

m - 3 (Feasible region is the unbounded region 
below the indicated nonlinear constraint manifold) 


Finally, the algorithm has two stopping conditions. 
the search is stopped if the change in the cost function and the 
change in the length of the independent-variable vector between 
two successive iterations fall below their respective inpu 
tolerances; that is, if 


' F (vi) - F M 


(VIII -7 5) 


•Vi - V ' E2 * 

c« ffMV < t k # procedure is discontinued if the number of the cur- 
rent Iteration equals the maximum permissible number input by 
the user. 
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